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This  report  investigates  the  scattering  of  a  monochromatic  acoustic  wave 
by  an  assumed  unchanging  atmospheric  turbulence  distribution  that  moves 
with  a  uniform  constant  horizontal  wind.  The  source  and  detector  are  at 
rest  at  different  positions  in  a  ground-fixed  frame.  Theoretical  expressions 
are  derived  for  the  scattered  acoustic  pressure  at  the  detector.  Because  dis¬ 
tances  and  scattering  angles  change  with  time  as  any  turbulent  fluctua¬ 
tion  (turbule)  flows  through  the  scattering  volume,  the  detector  signal  has 
a  time-dependent  Doppler  frequency  shift.  A  simple  formula  for  the  con¬ 
sequent  frequency  broadening  is  derived.  A  computer  code  is  developed 
that  calculates  the  scattered  signal  and  its  fast  Fourier  transform  scattered 
by  a  single  turbule,  or  by  a  collection  of  turbules  of  many  different  scale 
lengths.  The  code  uses  a  new  approximate  but  quite  accurate  time-shift  al¬ 
gorithm  that  reduces  the  calculation  time  substantially.  Several  numerical 
results  from  this  code  are  presented  for  a  model  scenario  that  is  similar  to  a 
recent  experimental  scenario,  with  good  agreement  among  spectral  widths 
and  shapes. 
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1.  Introduction 


This  report  investigates  acoustic  scattering  by  atmospheric  turbulence  that 
is  totally  frozen  (stationary)  in  a  frame  comoving  with  a  uniform  constant 
horizontal  wind.  The  ambient  temperature  is  also  assumed  constant  and 
uniform.  For  convenience,  the  turbulence  is  modeled  as  an  appropriate 
collection  of  self-similar  eddies  or  "turbules,"  as  described  in  past  litera¬ 
ture  [1].  A  turbule  is  a  localized  region  of  temperature  and/ or  flow  velocity 
fluctuations  with  respect  to  ambient  values. 

An  acoustic  detector  and  a  monochromatic  source  are  at  rest  in  a  ground- 
fixed  frame  at  different  positions.  The  detector  and  source  are  assumed  to 
be  in  the  far-field  region  of  a  turbule  at  any  instant  but  in  the  near-field 
region  of  the  complete  stationary  volume  that  the  turbules  traverse.  This 
is  almost  always  the  case  in  experimental  situations.  The  theory  developed 
during  this  work  first  obtains  a  Born  approximation  expression  for  the  time 
dependence  of  the  acoustic  pressure  at  the  detector  caused  by  scattering 
of  the  source  pressure  wave  by  one  turbule.  A  computer  algorithm  devel¬ 
oped  during  this  work  calculates  this  expression.  This  algorithm  also  calcu¬ 
lates  the  fast  Fourier  transform  (FFT)  of  this  pressure  wave  and  reveals  that 
the  wave's  frequency  spectrum  is  both  shifted  and  broadened  appreciably 
around  the  source  frequency.  The  broadening  occurs  because  the  scattering 
angles  and  distances  for  any  turbule  are  slowly  varying  functions  of  time 
for  realistic  wind  speeds.  A  simple  formula  for  the  received  frequency  as 
a  function  of  time  is  also  derived  by  two  different  methods;  it  predicts  the 
same  broadening  as  the  computer  algorithm. 

The  conventional  theory  [2]  does  not  allow  for  this  kind  of  time  depen¬ 
dence,  because  it  assumes  that  the  detector  and  source  are  in  the  far-field 
region  of  the  whole  scattering  volume.  Thus  it  predicts  spectral  broadening 
only  if  there  are  time-dependent  fluctuations  in  the  comoving  frame.  Estab¬ 
lishing  the  amount  of  spectral  broadening  and  shifting  for  turbulence  com¬ 
pletely  frozen  in  the  comoving  frame  but  for  realistic  near-field  detector 
and  source  locations  thus  became  one  of  the  most  important  objectives  and 
results  of  this  work.  The  experimentally  observed  spectral  broadening  used 
for  comparison  has  a  fairly  narrow  large  peak  around  the  source  frequency 
that  is  very  jagged  but  overall  appears  Gaussian,  with  wide  noisy-looking 
wings  at  smaller  amplitude.  The  theory  developed  in  this  report  predicts 
that  the  large  narrow  peak  is  due  to  large  turbules,  but  the  wide  wings  are 
due  to  small  ones,  and  also  predicts  that  the  jaggedness  is  due  to  random¬ 
ness  in  turbule  locations  in  the  comoving  frame.  The  computer  algorithm 
based  on  the  theory  yields  spectra  that  exhibit  the  same  general  features  as 
some  experimental  data,  including  the  spectral  widths  and  jaggedness,  for 
realistic  populations  of  turbules  of  many  length  scales. 
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Also  developed  in  this  work  is  an  efficient  approximate  method  for  adding 
contributions  of  many  small  turbules  of  each  size  that  are  approximately 
in-line  with  the  wind.  The  method  is  demonstrated  to  be  accurate,  and 
it  reduces  computer  running  time  significantly,  enough  so  that  three- 
dimensional  modeling  is  indeed  practical. 


2.  Technical  Discussion 


2.1  Scattering  by  One  Turbule 

2.1.1  Geometry 

The  scattering  geometry  is  defined  by  figure  1.  The  letters  S.  D.  and  T  stand 
for  source,  detector,  and  turbule,  respectively.  Coordinate  2  is  vertical,  z  =  0 
is  the  ground  plane,  and  coordinates  x  and  y  are  horizontal,  ground-fixed. 
The  uniform  constant  horizontal  wind  is  w,  given  by 

w  =  w(ex  cos  <pw  +  ey  sin  4>%v)  ■  (1) 

where  (ex.ey.ez)  are  the  Cartesian  unit  basis  vectors  and  <pu,  is  the  az¬ 
imuthal  angle  of  the  wind  as  shown.  The  source  S  and  detector  D  are  at 
locations  given  by  position  vectors  (R5,  Rc),  respectively  (fixed  with  re¬ 
spect  to  the  ground-fixed  reference  frame)  and  the  turbule  center  T  has 


time- varying  position  vector  Rr(t). 

The  location  of  the  "center"  of  the  turbule  T  is  given  by 

R T(t)  =  b  +  wt .  (2) 

so  b  is  the  location  of  the  turbule  at  t  =  0  in  the  ground-fixed  frame. 

The  vectors  R 7  ,s',  R ;  1 ;  are  defined  by 

Rrs(t)  =  R T(t)  -  R5,  R DT(t)  =  RD  -  R T(t) .  (3) 

In  the  comoving  (barred)  frame,  defined  by 

x  =  x-  wxt,  y  =  y-  wyt,  z  =  z  ,  (4) 

so  that  the  origins  coincide  at  t  =  0,  we  have 

RT  =  b,  Rs(f)  =  Rs  -  wt ,  RD{t)  =  RD-wt  (5) 

so  that 

R TS(t)  —  Rt  -  R s(t)  =  b  —  (RS  -  w t)  =  R TS{t)  and  (6) 
R DT{t)  =  R°{t)  -  Rr  =  Rd  -  wt  -  b  =  RDT(t)  .  (7) 


Figure  1.  Geometry  in 
ground-fixed  reference 
frame. 
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2.1.2  Scattered  Acoustic  Pressure  at  Detector 


We  seek  expressions  for  the  acoustic  pressure  waves  plas(t))  scat¬ 

tered  by  the  turbulent  (temperature,  velocity)  fields  of  the  turbule,  as  func¬ 
tions  of  time  at  the  detector.  These  expressions  are  derived  in  detail  in  ap¬ 
pendix  A.  Here,  we  quote  the  results  valid  to  all  orders  of  ,1  =  ir/r,  where  c 
is  the  ambient  adiabatic  sound  speed,  taken  to  be  constant  and  uniform  in 
the  scattering  volume  in  this  work.  The  source  is  assumed  to  be  monochro¬ 
matic  and  isotropic  in  its  rest  (ground-fixed)  frame.  Reflections  from  the 
ground  are  ignored.  From  equation  (A-71),  we  have  for  pj,'s(t)  the  form 

p!,M)  =  -4 '  (0  exp (-/u;,,/)  +  c.c.  .  (8) 

where  uos  is  the  source  angular  frequency  and  A1  (/)  is  a  (slowly  varying) 
complex-valued  amplitude,  given  by  equation  (A-72): 


AT(t)  =  -V, 


■RDT)T(K)exp[ikA  +  Hnr)]  . 


(9) 


where  a  hat  over  a  vector  indicates  that  it  is  a  unit  vector.  Here,  ks  =  uss/c, 
V,  is  the  pressure  of  the  source  wave  a  distance  a1  from  the  source  in  a  non¬ 
moving  atmosphere,  T0  is  the  constant  uniform  ambient  temperature,  and 

T(K)  is  the  spatial  FFT  of  T(r  —  b),  the  stationary  temperature  variation  of 
the  turbule  in  the  comoving  frame,  given  by 

f(K)  =  jd^TiOeM-iK-O-  (10) 


Also,  the  various  vectors  and  distances  are  defined  by  equations  (A-65)  to 
(A-70),  as  follows: 


where 


with 


*  M,,w  '"}• 

(11) 

as"  =  W’Vv  +  ■ 

(12) 

II 

vg 

1 

(13) 

i?;  =  njs(t)  - 

(14) 

£>*  =  [(i  -d-)m2-(p-RN)‘2Y2 

(15) 

4  =  =  y  {Ds  —  /3  •  R^} 

(16) 

T,  =  =  (l-d2)-^ 

(17) 

and  with  (R TS(t),  KPI  ( t ))  given  by  equations  (6)  and  (7).  Note  that  all  the 
quantities  in  equations  (10)  to  (16)  and  thus  AT(t)  are  functions  of  time  that 
vary  slowly  compared  to  exp (—iajJ). 
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These  results  are  quite  complicated,  primarily  because  of  retardation  ef¬ 
fects.  Only  when  /3  — *•  0  are  they  simple.  Then  we  get 

R0j  a  RSJ  «  Rjs.  Ds  »  Rs  ps  i?T5.  (18) 

A',  S3  k8(Ri>T  -  RJ's)  .  and  (19) 

At  =  -P,(^O74>rT„)(^^rs)f(K)exi.[a',(fl,'s  +  RDT)].(20) 

where  6  is  the  scattering  angle,  given  here  by 

cos  0  =  Rdt  ■  Rts  -  (21) 


These  quantities  are  still  weakly  time  dependent  for  d  <<<  1.  If  we  define 


ks  =  Rrs  ,  Apu,  =  -Vea!(RTS)  1  exp (ik.RTS)  . 


r  =  Rdt 


(22) 


then  the  amplitude  .4 7  of  equation  (20)  has  the  form 


Ar  =  Arm(?S^l) 


r 


(23) 


with 


cos  9  =  ks  •  r  ,  K  =  ks  (r  —  ks)  . 


(24) 


This  is  exactly  the  form  of  the  Monin  [3]  result  for  temperature  scatter¬ 
ing  from  stationary  turbulence,  where  Apw  is  the  amplitude  of  an  incident 
plane  wave  traveling  in  direction  ks. 

From  equation  (A-81),  the  contribution  p£a(f)  scattered  by  the  solenoidal 
turbulent  velocity  field  of  the  turbule  is 

Pas(t)  =  ^(f)  exp(-iust)  +  c.c. ,  (25) 

where  the  complex  amplitude  Av(t)  is  defined  by 


Av(t)  = 


kgVea'\  f  R^Rq  •  Rgr)(Rpr-  v  (K)) 
2ttc72  )  \  ( DS)A  RDT 


exp[ifcs(i?o  T  A/7^)]  . 


(26) 


with  v  (K)  the  spatial  FFT  of  the  turbule  velocity  field  v(r  —  b),  stationary 
in  the  comoving  (barred)  frame,  given  by 

v  (K)  =  Jds£exp(-iK  ■  ^)v(^)  .  (27) 

For  0  — >  0,  Av(t)  also  reduces  to  the  Monin  result  for  scattering  by  turbulent 
velocity  fluctuations. 

Equations  (9)  and  (26),  with  the  accompanying  equations  (10)  to  (17)  and 
(27),  are  used  in  the  computer  algorithm  described  in  appendix  B.  To  ob¬ 
tain  T(K)  and  v  (K)  specifically,  we  must  use  a  specific  turbule  model 
described  in  the  following  section. 
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2.1.3 


Turbule  Model 


For  convenience,  we  use  Gaussian  envelopes  and  spherical  symmetry,  as  in 
a  past  paper  [1],  A  temperature  turbule  of  scale  length  or  "size"  a,  fixed  in 
the  comoving  (barred)  frame,  is  written 


7'(r  —  b)  =  AT),  exp 


(f  —  b)2 


ci- 


and  a  velocity  turbule  of  size  «  is  written 

v(r  -  b)  =  [n„  x  (r  -  6)]  exp 


(r~b)J 


(28) 


(29) 


Note  that  ftn  is  an  angular  velocity,  and  v  is  solenoidal. 


The  amplitudes  (A Tn.  fif,  =  |fi„|)  are  scaled  with  a,  as  in  a  previous  paper 
[1].  Specifically,  we  designate  a l  as  the  size  of  the  largest  turbule  to  be  used, 
the  outer  scale  length  for  this  problem,  and  choose  values  for 


in  C 


(30) 


Typical  values  are  rr  «  T,.  «  0.01.  Here,  ST,  =  11,,,,  AT;  =  A T„t .  (We  must 
also  choose  the  direction  fla.)  Then,  from  the  previous  paper  [1]  and  also 
in  accordance  with  conventional  models  [2],  we  choose  for  any  size  a 


A Tn  =  (a/a^AT,  .  Slaa  =  (o/aj^.a,  ,  (31) 


or 


na  =  (aja)2^  . 


(32) 


From  the  definitions  of  equations  (10)  and  (27)  and  from  equations  (28)  and 
(29),  we  get  the  same  expressions  as  in  appendix  A, 


T{ K)  =  7T3/2  o3  AT),  exp  — - 

4 

v  (K)  =  —7 r3/,“o'r>ST0  x  Kexp 


(33) 

(34) 


Note  that  there  is  no  velocity  scattering  if  Qn  ||  K  . 


Consider  the  lowest-order  approximation  for  the  scattering  vector  K,  i.e., 
just  like  equation  (19): 

K (/)  «  ks(R DT{t)  -  Rrs(0)  •  (35) 


This  is  a  good  approximation  if  3  <§C  1.  Then 

A'2  =  2A'2(1  -  cos  (9(0)  =  4A-2  sin2  .  (36) 

with  the  scattering  angle  0(f)  defined  by  equation  (21).  Therefore,  if  ksa  1, 

then  T(K)  and  v  (K)  are  extremely  small  unless  6  is  near  0°,  so  the  scatter¬ 
ing  is  negligible  at  large  angles.  But  for  ksn  «  1,  the  Gaussian  envelope  in 
equations  (33)  and  (34)  is  essentially  unity  for  all  scattering  angles,  so  back¬ 
ward  scattering  is  comparable  to  near  forward  scattering  in  magnitude. 
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2.1.4  Spectral  Broadening 

We  note  that  both  AT(t)  and  Av(t)  contain  the  factor 

exp{— iujs[t  -  ( Ro(t )  +  RDT(t))/c}}  =  exp (-*$(f))  .  (37) 

where  the  equality  defines  the  phase  <J>.  Consider  any  time  t  =  t„.  Then 
write  for  t  near  fn 

$(*)  »  $(f0)  +  ($(f0))(*  -  *0)  +  neglect .  (38) 

where  we  can  neglect  terms  of  higher  order  in  (t  —  tu)  if  Ro{t)  and  Rm  it) 
are  "slowly  varying,"  which  is  the  case  for  0  <C  1.  Thus  4>(f0)  is  the  effective 
angular  frequency  u>(t0)  at  time  t0.  Since  t0  is  any  time,  we  have 

0j(t)  =LUS{  1  -  -c(R0(t)  +  RDT(t0)))  .  (39) 

When  equations  (12),  (3),  and  (13)  to  (17)  are  used,  obtaining  an  expression 
for  io(t)  valid  to  first  order  in  0  is  very  simple.  Since  c-1  times  the  time 
derivative  of  any  position  vector  or  distance  is  proportional  to  0,  then  for 
equations  (13)  to  (17),  we  may  put,  to  order  0, 

R0(t )  «  RTS{t )  =  dt{ Rts  ■  RTS)1/2  =  RTS  ■  RTS  =  w  •  Rts  (40) 
RDT(t)  =  dt(RDT  ■  Rdt)R'2  =  Rdt  ■  Rdt  =  -w  •  Rdt  .  (41) 

Therefore,  to  order  0, 

u(t)  «  US[1  +  0  ■  (R DT{t)  -  R TSm  .  (42) 

Note  that  this  can  be  written  in  terms  of  the  approximate  scattering  vector 
K  (t)  of  equation  (35)  as 

u(t)  =  u>s  +  w  •  K(t)  .  (43) 

This  relation  is  actually  correct  for  the  exact  K (t)  given  by  equation  (11). 

From  equations  (32)  or  (33),  (36),  and  (33)  and  (34),  and  from  figures  1(a) 
and  1(b),  it  is  evident  that  small  turbules  may  produce  much  wider  spectral 
broadening  than  large  ones.  This  is  because  the  Gaussian  envelope  cuts  off 
the  scattering  from  large  turbules  when  |w  •  K|  is  "large,"  but  does  not 
do  so  for  small  turbules.  Thus  from  small  turbules,  we  can  get  appreciable 
backscattering,  so  the  relative  frequency  shift  w  •  Kws  may  vary  from  zero 
to  ±20,  while  for  large  turbules,  |w  •  K|ws  <  20  for  any  scattering  angles 
that  produce  appreciable  scattering. 

It  should  be  emphasized  that  this  general  behavior  is  not  peculiar  to  a  tur- 
bule  model,  or  to  Gaussian  model  turbules,  but  is  common  for  any  model 
of  turbulence  that  contains  both  large  and  small  length  scales  and  is  frozen 
in  a  comoving  frame,  provided  only  that  the  source  and  detector  are  in  the 
near-field  region  of  the  complete  scattering  volume  Vs.  This  applies  to  the 
scenarios  of  interest  in  this  work.  If  (and  only  if)  the  source  and  detector 
are  both  in  the  far-field  region  of  Vs  at  all  times,  then  the  conventional  re¬ 
sult  ensues  that  moving  frozen  turbulence  may  produce  only  a  frequency 
shift  but  no  broadening  [2]. 

In  appendix  C,  a  different  simple  derivation  of  equation  (42)  is  provided. 
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2.2  Multiturbule  Model 


This  section  contains  a  summary  description  of  the  scheme  used  later  in 
the  report  for  populating  the  scattering  volume  with  turbules.  Rather  than 
having  a  continuous  distribution  of  sizes  for  the  turbule  size  parameter  a,  a 
discrete  set  of  sizes  «„  will  be  selected.  For  each  size  class  index  v ,  a  cubical 
cell  containing  one  turbule  is  defined  whose  linear  dimension  d„  is  an  ap¬ 
propriate  multiple  of  nn.  Suppose  that  the  time  interval  for  which  data  are 
desired  is  T  and  that  the  average  wind  velocity  is  w.  Then  an  approximate 
scattering  volume  length  L  is  equal  to  wT.  A  tube  is  defined  to  be  a  re¬ 
gion  of  space  consisting  of  a  sufficient  number  of  cells  AI„  aligned  along  the 
wind  direction  so  that  the  total  length  of  the  tube  is  L  or  greater.  The  several 
tubes  for  each  size  class  n  are  stacked  laterally  and  vertically  to  overspread 
the  scattering  volume.  For  this  discussion,  consider  that  a„+\/a„  —  1/2. 
Then  dn+\/d„  =  1/2  and  Mn+\/M„  =  2.  The  number  of  tubes  for  class 
■n  +  1  is  four  times  the  number  for  class  v. 

2.2.1  General  Population 

Our  interest  is  in  the  scattered  acoustic  pressure  produced  by  a  multitur¬ 
bule  model  that  contains  appropriate  numbers  of  turbules  of  many  sizes, 
from  the  largest  (outer  scale)  a  t  to  the  smallest  (inner  scale)  aA,  ,  with  Ns 
different  scales  included.  As  shown  by  Goedecke  et  al  [1],  if  we  define  the 
packing  fraction 


(t>  -  n„o 3  .  (44) 

where  na  is  the  number  density  of  turbules  of  size  a,  then  </>  must  be  inde¬ 
pendent  of  a  so  that  the  energy  transfer  rate  per  unit  mass  in  a  Kolmogorov 
cascade  can  be  scale-independent  as  well.  On  average,  one  turbule  of  size 
a  is  somewhere  in  each  cubical  cell  of  size  d„,  where 

d«  =  71  a  1/3  =  .  (45) 

For  a  typical  <j>  %  10~3  ,  da  =  10a  . 

If  we  had  unlimited  computer  resources,  we  could  populate  a  defined  scat¬ 
tering  volume  Vs  with  turbules  of  different  sizes  (one  "randomly"  located 
in  each  cell  for  each  size)  and  simply  add  up  the  ( AT (t ).  Al’(t))  from  each. 
This  is  actually  what  we  do  for  the  larger  turbules  with  ksa  >  1.  But  if  there 
are  perhaps  100  turbules  with  a  =  2  m  in  VS/  then  for  a  =  0.02  m  =  2  cm,  ac¬ 
cording  to  equation  (44),  there  will  be  (100)(2/0.02)3  =  108  turbules,  a  truly 
formidable  number. 

2.2.2  Time-Shift  Algorithm 

Fortunately,  there  seems  to  be  a  good  approximate  way  to  handle  the  mul¬ 
titude  of  smaller  turbules  that  reduces  the  calculation  time  significantly. 
Imagine  two  identical  turbules,  1  and  2,  more  or  less  in-line  along  the  wind 
direction,  chosen  to  be  the  ./  -direction.  That  is,  let  the  centers  of  their  two 
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cells  have  coordinates  (b&,  bq),  b^z)  and  (b&,  bqj ,  b^J);  i.e.,  their  cells  are 
in-line  with  the  wind. 

Now,  consider  equation  (9)  for  from  turbule  1,  and  assume  that  the 

source  and  detector  are  always  in  the  radiation  zone  of  a  comoving  cell  of 
side  d  =  in  which  the  turbule  is  located.  This  implies  R^J (t)  and 

Rji){f)  >  (d.  A.  d2/ A),  whichever  is  greater,  at  all  (or  most)  times.  For  such 

cases,  vectors  such  as  (Rr5^\  . . . )  can  be  replaced  by  (Rj5^. . . . ),  the  vec¬ 
tors  to  the  center  of  the  cell,  except  in  the  factor  exp [iks(R^  +  BDT<i))}.  In 
this  factor,  we  may  write 

exp [iks{R[V  +  Rdt{1))}  «  exp[iks(R^  +  i?fT(1))]  exp(-iK^1)  •  r^)  ,  (46) 

where  77  x  is  the  displacement  of  the  center  of  turbule  1  from  its  cell  center. 
So  we  get 

AJd (t)  =  exp(— iK^(f)  •  Vi)  ,  (47) 

where  from  equations  (9)  and  (46),  we  have 

Aj^t)  =  -Pe(fefa/47rToo72)[(41c))2/(4S(1))3^r(1)] 

x(il1c)-  KfT^)T0(K^)eMiks(R{oJ  +  R?1'™)]-  (48) 

That  is,  Aj^  (t)  is  the  old  AjV)  (t)  but  calculated  with  vectors  to  the  cell  center 
times  the  phase  factor  exp[iks(Rt0]}  +  R®Ti1^)]. 

Now,  for  turbule  2,  clearly  (Aj2^(t),  Ki2)(f))  are  simply  time-shifted  from 
(*4fi)(*)’  K^(f)).  That  is,  we  have 

A[2)(t)  =  AT(1){t  +  r2) ,  K f  (*)  =  +  r2)  ,  (49) 

where  r2  =  w~l(b^  -  b&).  For  example,  suppose  turbule  2  is  one  cell 
ahead  of  turbule  1,  so  that  b^J  —  b^J  =  d.  Then  r2  =  d/w,  so  Aj2){t )  is  just 
like  Aj^ (t),  but  occurs  earlier. 

Thus,  the  time-shift  algorithm  for  a  whole  tube  of  M  turbules  of  a  given 
size,  with  cell  centers  equally  spaced  by  d  along  the  wind  direction,  is 
obviously 

AI 

^tabe(*)  =  Ajx)  (t  +  T()  exp[-iK<.1}  (t  +  T()  ■  Ve\ ,  (50) 

e=i 

where 


n  =  w-'(blg-V»)  =  (e-l)(d/w),  (51) 

and  T](  is  chosen  "randomly"  for  each  £;  i.e.,  —d/4  <  {iiexAltyVCz)  <  d/4, 
with  uniform  independent  probabilities.  So,  for  a  given  tube,  we  need  to 
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compute  the  complicated  functions  Kc'^/))  only  once;  then,  for 

each  turbule  l  in  the  tube,  we  choose  an  r/(  and  compute 

Ajn(t  -  T( )  =  Ajt](t)exV(-Ki]\t)  ■  n,)  .  (52) 

then  simply  time-shift  this  to  get  A^(t),  and  add  that  to  the  accumulated 

By  comparing  with  direct  computation  of  the  A{uhe(t)  for  the  same  tube 
and  using  the  full  equations  (9)  to  (16)  for  each  turbule  with  its  center 
at  the  same  set  of  random  locations  as  rjf  from  above,  we  found  that  for 
ia,  =  500  Hz,  the  time-shift  algorithm  is  a  good  approximation  for  d  < 
4  m.  However,  the  algorithm  is  not  as  good  for  d  >  8  m,  with  R?r  and 
Rjs  >  10  m  minimum  but  actually  >  100  m  during  the  times  of  largest  sig¬ 
nal  production.  Since  a  direct  computation  involves  ~100  operations  per 
turbule  and  a  time-shift  of  only  ~5  operations,  the  computing  time  advan¬ 
tage  of  the  time-shift  algorithm  is  major. 

These  same  considerations  apply  for  the  Av(t)  of  equation  (26),  except  that 
we  need  a  separate  ST  for  each  turbule.  For  any  turbule,  we  may  write  in 
general 


Av(t)  =  ft  ■  Dv(t)  ,  (53) 

where  from  equations  (26)  and  (34),  we  have 


Rl(K0-KDT)\ 
(. DS)4RDT  ) 


exp(— A'2o2/4)  exp[iks(Ri)  +  RDT)](RD1  x  K)  . 


(54) 


The  time-shift  algorithm  for  a  tube  of  M  turbules  along  the  wind  is  now 

m 

Kxbei*)  =  K)  ^t  +  T^  pxp[-':kc  }  (f  +  n )  ■  V(]  ■  (55) 

(=\ 

where  (B^\  K*1 )(/))  are  calculated  with  vectors  to  the  center  of  the  cell 
containing  and  comoving  with  turbule  1.  We  calculate  the  same  way  as  for 
A[uhe(t),  except  that  we  must  choose  an  S %  as  well  as  an  r/(  for  each  turbule. 


2.3  Results 


Many  of  the  theoretical  results  are  presented  and  discussed  in  the  previous 
section  and  in  appendixes  A  to  C,  so  we  limit  this  section  primarily  to  pre¬ 
senting  and  discussing  graphs  that  illustrate  some  of  the  most  interesting 
predictions  of  the  computer  code  (Slow  Eddy)  that  was  developed  on  the 
basis  of  the  theoretical  results.  However,  we  will  also  show  how  the  theory 
easily  predicts  the  major  features  of  these  computational  results. 
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For  all  the  following  results,  the  detector  was  placed  at  (—150,  0,  0)  m  and 
the  source  at  (150,  0,  10)  m,  and  the  wind  was  chosen  with  components 
(3.1,  0,  0)  m/s.  Thus,  the  primed  and  unprimed  ground-fixed  frames  de¬ 
scribed  in  section  2.1.1  and  discussed  in  appendix  B  were  the  same,  and 
the  wind  velocity  was  blowing  from  the  detector  toward  the  source,  in  the 
(+,r)  direction.  The  total  time  T  was  24  s,  0  <  t  <  T,  tube  length  was 
L  =  v’T  =  756.4  m  centered  at  the  origin,  the  time  steps  At  were  0.0403  s, 
and  the  total  number  of  points  used  in  the  calculation  was  N  —  6048.  (The 
point  n  =  6049  <-+  t  =  T  is  "equivalent"  to  t  =  0  because  of  the  inher¬ 
ent  periodicity  of  the  discrete  FFT.)  The  adiabatic  sound  speed  was  taken 
to  be  c  =  340  m/s,  corresponding  to  an  ambient  temperature  T0  «  15  °C 
«  288  K  «  60  °F.  For  multiturbule  computations,  each  turbule  angular 
velocity  unit  vector  was  generally  chosen  to  have  uniform  random  distri¬ 
bution  in  47T  solid  angle,  independently  of  other  turbules.  This  orientation 
distribution  produces  a  model  of  isotropic  turbulence,  as  was  discussed  in 
another  paper  [1].  The  A Ta  in  the  scattering  volume  Vs  were  generally  cho¬ 
sen  to  have  the  same  sign  (positive)  because  temperature  turbules  are  air 
pockets  warmer  than  ambient  that  have  risen  from  the  ground  because  of 
buoyancy  forces,  such  as  those  over  a  paved  runway  that  was  used  in  the 
relevant  experiment.  Over  a  larger  volume,  the  mean  AT  probably  should 
be  zero.  In  the  following  results,  one  calculation  was  made  with  each  tur¬ 
bule  AT  either  >  0  or  <  0  with  equal  probability  to  contrast  the  results  of 
the  different  conditions.  Also,  some  calculations  were  made  with  the  an¬ 
gular  velocity  vectors  of  the  larger  turbules  all  in  the  -by-direction,  for  the 
same  reason. 

The  ratio  of  cell  size  to  turbule  size  was  chosen  to  be  d/a  —  8,  yielding 
a  packing  fraction  <j>  =  1.95  x  10-3.  The  number  of  sizes  in  the  cascade 
was  chosen  to  be  Ns  =  6,  and  the  (largest,  smallest)  turbule  sizes  (a15  aNs) 
were  chosen  as  (1.5  m,  0.046825  m),  respectively.  Thus  the  parameter  Nd 
of  equation  (B-18)  was  equal  to  2  (whereby  q  =  In  2),  such  that  each  suc¬ 
cessively  smaller  size  was  half  the  previous  size.  The  relative  fluctuation 
strengths  (ATi/T0  =  Tt.  VLlal/c  =  I\,)  of  equation  (30)  were  chosen  as 
(Tr  =  0.013,  I\,  =  0.027).  From  our  previous  paper  [1]  or  equation  (B-17) 
of  this  report,  these  values  yield  C2  =  0.088  m4/3/  s2.  C2  =  0.044  K2/m2'/3. 
These  values  are  in  good  agreement  with  those  given  by  Brown  and  Clif¬ 
ford  [4].  Also,  this  value  of  C'2  was  inferred*  from  wind  data  that  were 
measured,  along  with  the  time-dependent  acoustic  pressure  received  by 
many  different  detectors,  during  an  experiment  by  Galindo  and  Havelock 
[5].  All  the  numerical  values  chosen  in  this  report  are  supposed  to  provide 
a  good  approximation  to  those  that  existed  in  that  experiment.  (As  men¬ 
tioned  earlier,  a  main  objective  of  this  work  is  to  determine  whether  a  real¬ 
istic  model  of  acoustic  scattering  by  turbulence  that  is  frozen  in  a  frame  and 
comoving  with  a  wind  can  possibly  simulate  actual  experimental  results.) 

Using  these  values  of  the  input  parameters  and  the  method  and  computer 
code  described  in  this  subsection  and  in  appendix  B,  we  calculated  total 

"D.  K.  Wilson  of  ARL  computed  this  value. 
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time-dependent  complex  acoustic  pressure  amplitudes  ( Ar(t ).  Av(t))  scat¬ 
tered  by  individual  turbules,  tubes  of  turbules,  and  all  tubes  of  turbules  in 
several  different  scattering  volumes  and  then  Fourier  transformed  to  get 
the  spectra  (A1  (a-').  Ar  (*>)).  For  convenience,  each  (Ar .  Av)  of  equations  (9) 
and  (26)  was  "normalized"  by  division  by  the  factor  V,  with  o'  —  1  m.  Thus, 
V,  is  now  the  pressure  amplitude  of  the  source  wave  at  a  distance  a'  =  1  nr 
from  the  source  in  a  still  atmosphere.  In  what  follows,  the  ( A 1 .  Ar)  are  the 
normalized  ones,  which  are  dimensionless. 

In  figure  2,  we  plot  the  real  and  imaginary  parts  and  the  magnitudes  of 
( AT(t ).  Ar(t)),  respectively,  for  a  single  turbule  with  a  =  o,  =  1.5  m, 
with  bj.  —  —L/ 2.  btl  =  0.  bz  —  20  nr.  Here,  b~  —  20  nr  corresponds  to 
an  above-ground  base  height  zi,  =  14  nr,  plus  r/,/2  =  6  nr,  half  the  side 
of  a  cubical  cell  having  <7,  =  8<j,.  In  figure  3,  we  plot  the  magnitudes  of 
the  FFTs  A1  {^)  and  .4r(^).  One  can  easily  understand  from  equation  (33) 
why  the  pulse  in  figure  2  has  a  fairly  narrow  width  in  time:  The  Gaus¬ 
sian  envelope  rapidly  kills  the  signal  as  the  scattering  angle  increases,  since 
kgO  j  ps  13.9.  The  frequency  width  can  be  obtained  easily  from  equation  (42), 
u>(t)  «  cj,[1  +  /3  •  (RD7  (/)  -  R75(/))] .  That  is,  calculate  or  obtain  from  figure 
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2  the  two  times  .  t2 )  at  which  the  time-dependent  magnitudes  |  A(t)  |  are 
perhaps  5  percent  oftheir  maximum  values  and  find  R/r/  and  RTS  at  these 
times.  Then  verify  from  figure  3  that  |  ui(ti)  —  w(t 2)  |  is  a  good  estimate  of 
the  (angular)  frequency  width  of  the  A(u).  This  calculation  works  because 
the  complete  frequency  range  vs(  1  —  2d)  to  vs{  1  -I-  2/3)  is  "scanned"  by  any 
turbule,  but  only  a  small  part  of  it  is  manifested  in  the  scattered  waves  from 
large  turbules. 

In  figures  4  and  5,  we  repeat  the  presentations  of  figures  2  and  3,  respec¬ 
tively,  but  for  a  single  small  turbule  of  size  a6  =  4.6875  cm,  with  the  same 
(i bx .  by.bz).  For  this  size  turbule,  ( ksaG )2  =  0.188,  so  the  Gaussian  envelope 
of  the  A(t)  is  almost  unity  for  all  scattering  angles.  Therefore,  each  .4(f)  is 
nonzero  over  the  total  time  duration,  and  the  frequency  width  extends  to 
almost  the  maximum  possible,  ±2 0u3,  obtainable  by  backscattering  from 
equation  (42).  Note  also  that  each  A(u)  goes  to  zero  at  frequencies  = 
4.3  Hz,  v-  =  -4.2  Hz.  Again,  this  is  easily  understood  from  equations  (9), 
(26),  and  (42):  The  time-dependent  amplitude  A(t.)  is  zero  approximately 
when  Rdt  ■  RTS  =  cos  0  =0  (to  lowest  order  in  0).  From  figure  3,  ob¬ 
tain  the  times  t.3,  t4  at  which  this  is  true,  and  then  use  equation  (42)  to  get 
and  u)(t4);  these  agree  quite  well  (within  5-10%)  with  the  frequencies 
v±.  (The  algebraic  results  are  very  sensitive  to  (f3,  f4)  values,  and  the  com¬ 
puter  algorithm  contains  higher  order  terms  in  0  that  slightly  influence  the 
results.)  Also  note  the  asymmetry  of  the  A(t)  and  A(u>);  this  is  due  to  the 
source  height  being  different  from  the  detector  height,  in  this  geometry. 

Figures  6  to  9  contrast  the  effects  of  randomness  versus  regularity  in  turbule 
center  locations  for  many  turbules.  In  figure  6,  we  present  |  Av(t)  |  for  a  tube 
of  turbules  with  a  =  0.75  m,  with  all  angular  velocities  in  the  y-direction 
and  each  turbule  center  at  the  center  of  its  comoving  cell  of  side  length 
da  =  8a  =  6  m.  This  tube  forms  a  moving  "diffraction  grating."  The  tube 
centerline  is  ( by  =  0,  bz  =  20  m),  and  there  are  Ma  =  L/d,a  =  126  turbules 
contributing  at  any  time.  (In  all  calculations  done  for  tubes,  each  cell  in  a 
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Figure  4.  Pressure  signal  (a)  3x10-9 
amplitude  scattered  by 
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fixed  tube  of  total  length  L  is  kept  filled;  i.e.,  as  a  turbule  is  leaving  at  L/2, 
another  is  entering  at  -L/2.  This  ensures  steady-state  conditions.)  The  time 
dependence  is  extremely  regular,  as  one  would  expect. 

Figure  7  is  for  the  same  conditions  as  figure  6,  except  that  each  turbule 
center  is  randomly  placed  inside  its  comoving  cell,  with  no  correlations 
among  placements  of  different  turbules.  The  randomness  consists  of  uni¬ 
form  uncorrelated  probabilities  of  a  turbule  center's  (x.  y.  z )  locations  being 
between  -d/A  and  d/A  from  its  comoving  cell  center.  It  is  interesting  (but 
totally  expected)  that  this  randomness  produces  extremely  large  relative 
fluctuations. 
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Figure  8.  Signal 
frequency  spectrum 
scattered  by  velocity 
fluctuations  for  one 
tube  of  0.75-m  turbules, 
equally  spaced,  with 
spins  in  same  direction. 


Figure  9.  Signal 
frequency  spectra 
scattered  by  velocity 
fluctuations  for  one 
tube  of  0.75-m  turbules, 
randomly  spaced,  with 
spins  in  same  direction. 


In  figures  8  and  9,  we  plot  the  magnitudes  of  the  FFTs  for  the  conditions  of 
figures  6  and  7,  respectively,  on  a  semilog  plot  for  good  visibility.  In  figure 
8,  the  principal  maxima  are  striking;  they  are  spaced  just  as  they  should 
be,  i.e.,  at  frequencies  vq  =  qMa/T  «  (q)  (0.510)  FIz,  where  q  is  an  inte¬ 
ger.  There  should  be  exactly  zero  contribution  for  other  frequencies,  for 
the  discrete  FFT.  The  principal  maxima  diminish  rapidly  as  \v\  increases, 
probably  because  of  the  Gaussian  envelope  of  a  single  turbule.  (The  pic¬ 
ture  is  analogous  to  a  multiple-slit  interference  pattern  as  modified  by  the 
single-slit  diffraction  pattern,  in  optics.)  In  figure  9,  the  randomness  in  tur¬ 
bule  placement  yields  a  jagged  continuous  spectrum.  Both  figures  8  and 
9  display  a  remarkable  straight-line  behavior,  which  corresponds  on  these 
semilog  plots  to  an  exponential  spectral  envelope,  exp(-|^|r„),  where  ra  is 
a  constant  with  the  dimension  of  time  that  most  likely  depends  on  the  scale 
length  ci.  We  have  not  yet  been  able  to  predict  this  behavior  theoretically 
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Frequency  (Hz) 
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Figures  10  to  12  are  intended  to  illustrate  that  the  time-shift  algorithm  gives 
results  essentially  the  same  as  direct  calculation.  These  figures  are  \AT(u>)\ 
for  a  tube  of  Ma  =  504  turbules  with  a  —  1.5/23  =  0.1875  m,  with  the 
centerline  as  previously  mentioned,  with  all  A Ta  >  0.  Figure  10  is  for  reg¬ 
ular  placement  (for  reference),  with  the  use  of  "either"  the  direct  or  the 
time-shift  calculation;  i.e.,  the  results  are  identical  for  regular  placement 
of  the  turbule  centers  at  their  cell  centers  initially  Figure  11  is  for  random 
placements  using  the  direct  calculation.  Figure  12  is  for  the  "same"  random 
placements  but  using  the  time-shift  calculation,  which  is  eight  times  faster 
than  the  direct  method.  Clearly  the  two  methods  produce  substantially  the 
same  results,  although  there  are  small  differences,  as  there  should  be,  be¬ 
cause  the  time-shift  method  is  not  exact. 


Figure  10.  Signal 
frequency  spectrum 
scattered  by 
temperature 
fluctuations  for  one 
tube  of  0.1 9-m  turbules, 
equally  spaced,  using 
direct  method. 


0 

~o 

=3 

C 

Ui 

CO 

£ 


o 

0 

CL 

0 


D 

0 

0 

0 


> 

to 

0 

cc 


Figure  11.  Signal 
frequency  spectrum 
scattered  by 
temperature 
fluctuations  for  one 
tube  of  0.19-m  turbules, 
randomly  spaced,  using 
direct  method. 
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Figure  12.  Signal 
frequency  spectrum 
scattered  by 
temperature 
fluctuations  for  one 
tube  of  0.19-m  turbules, 
randomly  spaced,  using 
time-shift  method. 


Frequency  (Hz) 


Figures  13  and  14  illustrate  the  effects  of  random  versus  nonrandom  angu¬ 
lar  velocity  directions  fl  and  temperature  amplitude  signs  AT/|A7’|.  They 
are  for  a  tube  of  M„  —  63  turbules  with  a  =  1.5  m,  with  the  centerline  at 
by  =  6  m,  bz  =20  m.  Figure  13  plots  |^'  (u,')|;  figure  13(a)  is  for  all  the  ft 
in  the  //-direction,  while  figure  13(b)  is  for  the  ft  chosen  isotropically  at 
random.  Figure  14  plots  \A7\u;)\;  figure  14(a)  is  for  all  AT  >  0,  while  fig¬ 
ure  14(b)  is  for  sgn(AT)  =  ±1  with  equal  probability  (this  would  make 
T0  equal  to  the  mean  temperature  in  Vs).  In  each  case,  the  turbule  center 
locations  were  randomized  with  the  same  set  of  random  numbers,  so  dif¬ 
ferences  in  the  comparable  graphs  must  be  caused  by  the  different  choices 
of  ft  or  sgn(AT).  What  is  remarkable  is  that  the  differences  are  not  large. 
This  can  be  understood  by  noting  that,  for  example,  changing  the  sign  of 
AT/|AT|  for  a  given  turbule  introduces  a  phase  change  of  n  in  its  scat¬ 
tered  wave,  which  could  have  been  done  equally  well  by  moving  its  center 
about  a  half  wavelength.  So,  in  general,  randomization  of  ft  and  AT/\AT\ 
is  nearly  equivalent  to  nonrandomization  with  a  little  more  "freedom  of 
movement"  of  the  turbule  centers. 
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Figure  13.  Signal 
frequency  spectrum 
scattered  by  velocity 
fluctuations  for  one 
tube  of  1.5-m  turbules, 
randomly  spaced,  with 

(a)  spins  in  same 
direction  and 

(b)  random  spin 
directions. 
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Figures  15  to  17  are  sums  over  contributions  of  many  turbules  of  all  sizes 
used,  i.e.,  aa  =  (1.5,  0.75,  0.375,  0.1875,  0.09375,  0.046875)  m  for  a  -  (1, 
...,  Ns  —  6),  respectively.  The  amplitudes  (fla.  |AT„|)  were  scaled  as  dis¬ 
cussed  earlier,  so  these  calculations  are  an  attempt  to  numerically  represent 
realistic,  isotropic,  homogeneous,  fully  developed  turbulence  by  a  turbule 
model,  following  the  formulas  by  Goedecke  and  Auvermann  [1],  This  at¬ 
tempt  is  a  compromise  with  the  computing  power  of  a  PC  with  a  Pentium  II 
233  MHz  processor.  We  really  should  have  included  both  larger  and  smaller 
turbules,  plus  many  more  different  sizes  in  between.  But  what  we  did  al¬ 
ready  required  -4.25  hr  per  large  tube  (di  =  8«i  =  12  m),  corresponding 
to  30  s  per  turbule,  and  this  was  only  after  we  compromised  further  by 
making  a  reasonable  approximation  for  the  smallest  size  aG  =  0.046875  m. 
We  randomly  put  only  "one"  such  turbule  instead  of  eight  in  each  cell  of 
the  next  "larger"  size,  d5  =  8 a5  =  0.75  m,  and  then  we  multiplied  the  re¬ 
sult  by  eight.  We  checked  this  for  a  single  tube  of  side  length  d-,  versus  the 
"correct"  method  using  turbules  of  scale  aG  only,  and  the  results  were  quite 
close.  Since  so  many  such  tubes  are  in  one  big  tube  ((di/d5)2  =  256),  the 
total  error  induced  by  this  economy  should  be  quite  small. 

Figure  15  plots  |Ar(ic)  +  Au(w)|  for  two  big  tubes  with  centerlines  by  — 
±  6  m,  bz  —20  m.  Figure  16  does  the  same  for  two  big  tubes  above  those 
of  figure  15:  i.e.,  by  =  ±  6  m,  bz  =  32  m.  Figure  17  is  the  total  frequency 
spectrum  from  these  four  tubes.  Note  that  the  contributions  from  the  higher 
big  tubes  are  significantly  smaller,  and  that  the  general  shapes  of  the  plots 
are  essentially  unchanged,  so  these  shapes  are  established  already  by  just 
one  or  two  tubes.  The  semilog  plots  clearly  show  the  contributions  of  the 
small  turbules  all  the  way  to  ±2 /3vs  =  ±  10  Hz. 


Frequency  (Hz)  Frequency  (Hz) 

Figure  15.  (a)  Total  frequency  spectrum  due  to  temperature  and  velocity  fluctua¬ 
tions  for  two  full  tubes  with  sizes  1.5  m  down  to  4.7  cm.  (b)  Log  plot  of  (a)  from 
±10  Hz. 
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Figure  16.  (a)  Total  frequency  spectrum  due  to  temperature  and  velocity  fluctua¬ 
tions  for  two  full  tubes  located  above  the  two  tubes  in  figure  14,  with  sizes  1.5  m 
down  to  4.7  cm.  (b)  Log  plot  of  (a)  from  ±10  Flz. 
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Figure  17.  (a)  Total  frequency  spectrum  due  to  temperature  and  velocity  fluctua¬ 
tions  for  all  four  full  tubes  from  figure  14  and  15.  (b)  Log  plot  of  (a)  from  ±10  Hz. 
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Finally,  figure  18  contrasts  some  experimental  data  with  our  model  results. 
Figure  18(a)  is  from  a  presentation  by  Auvermann  and  Goedecke  [6];  the 
jagged  curve  is  the  data.  The  central  spike  is  about  101.56  Hz;  this  cor¬ 
responds  to  //  =  0  in  the  previous  spectral  plots  and  in  figure  18(b).  The 
ratio  of  the  central  spike  height  to  the  "mean"  height  at  v  =  102  -  101.56 
=  0.44  Hz  is  about  6/0.23  —  26.  These  data  were  for  a  source  detector  and 
wind  velocity  scenario  very  similar  (but  not  identical)  to  what  we  modeled. 

Figure  18(b)  is  an  expansion  of  our  figure  17(a)  along  the  abscissa  to  match 
the  total  frequency  width  in  figure  18(a).  The  model  predictions  are  simi¬ 
lar  to  the  experimental  results;  e.g,  the  ratio  of  the  central  spike  height  (at 
u  —  0)  to  the  "mean"  height  at  v  =  0.44  Hz  is  «0.52/0.02  =  26,  in  figure 
18(b).  But  the  jaggedness  in  the  model  results  is  not  as  great  as  in  the  exper¬ 
imental  results.  In  later  investigations,  it  has  been  found  that  the  jagged¬ 
ness,  defined  as  the  mean  frequency  spacing  of  adjacent  peaks  and  valleys 
in  the  spectrum,  is  approximately  equal  to  3 /T,  where  T  is  the  duration  of 
the  time  signal.  The  model  jaggedness  agrees  with  the  experimental  jagged¬ 
ness,  if  T  is  the  same  for  both. 


Figure  18.  Comparison 
of  experimental  data 
with  model:  (a)  Freq¬ 
uency  spectrum  from 
experiment  and  (b)  freq¬ 
uency  spectrum  from 
figure  17(a). 
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3.  Conclusions  and  Recommendations 


The  theory  developed  in  this  work  allows  understanding  of  the  spectral 
broadening  and  shifting  in  an  acoustic  wave  scattered  from  moving  tur¬ 
bulence  on  the  basis  of  two  simple  general  physical  principles.  The  first 
principle  is  that,  even  if  the  turbulence  is  frozen  in  a  frame  comoving  with 
a  wind,  spectral  broadening  occurs  because  the  scattering  angles  and  the 
distances  from  source  to  turbulent  eddies  and  from  the  eddies  to  the  de¬ 
tector  are  time-dependent.  The  analog  in  everyday  experience  is  what  you 
hear  when  an  approaching  car's  horn  is  blowing  at  frequency  us  and  then 
recedes  while  you  are  standing  to  one  side  of  the  car's  path:  You  actually 
hear  the  whole  Doppler-shifted  frequency  range,  from  //„.(  1  +  d)  to  us{l-li), 
as  a  function  of  time.  (For  scattered  waves,  the  factors  are  1  ±23,  where  /V  = 
(speed  of  object/ speed  of  sound),  to  first  order  in  d.) 

The  second  principle  is  that  large  objects  (compared  to  the  wavelength) 
scatter  strongly  but  mainly  in  the  near-forward  direction,  while  small  ob¬ 
jects  scatter  weakly  but  nearly  isotropically.  Thus  from  large  eddies,  we  get 
only  a  small  spectral  width,  while  from  small  eddies,  we  get  ±2d/vs  but  less 
amplitude. 

Of  course,  in  actuality  the  turbulent  flow  in  a  frame  comoving  with  the 
mean  wind  is  also  time-dependent.  But  we  have  shown  via  these  general 
principles  and  a  computer  model  that  this  time  dependence  is  unneces¬ 
sary  to  explain  the  general  features  of  some  observed  acoustic  data.  This 
is  markedly  in  contrast  to  the  conventional  treatment  of  spectral  broaden¬ 
ing,  which  attributes  all  of  spectral  broadening  to  time  dependence  in  the 
comoving  frame.  In  that  treatment,  the  source  and  detector  are  assumed  to 
be  in  the  far  field  of  the  whole  turbulent  scattering  volume,  so  that,  effec¬ 
tively,  scattering  angles  and  distances  are  not  time-dependent.  This  is  an 
unrealistic  assumption  for  acoustics,  where  both  the  source  and  detector 
are  invariably  in  the  near  field  of  the  scattering  volume. 

One  of  the  approximations  apparently  made  in  the  theory  developed  in  this 
report  was  the  assumption  of  a  "uniform"  mean  wind  and  a  "uniform"  am¬ 
bient  temperature.  (The  commonly  used  ambient  temperature  is  taken  as 
the  mean  temperature  in  the  scattering  volume  V«;  we  did  not  so  restrict  our 
ambient  temperature.)  But  these  may  only  be  "apparent"  approximations. 
For  example,  suppose  we  define  the  mean  wind  as  the  time  and  space  av¬ 
eraged  wind  over  the  time  of  an  experiment  and  over  the  whole  scattering 
volume  from  the  ground  up.  Then  of  course  it  is  by  definition  uniform  and 
constant,  so  "all"  wind  variations  are  part  of  the  turbulent  flow.  This  point 
of  view  would  introduce  anisotropy  and/or  inhomogeneity  in  the  turbu¬ 
lent  flow,  which  is  easy  to  do  in  a  turbule  model.  Also,  suppose  we  define 
the  mean  temperature  as  the  same  time  and  space  average  of  the  nona¬ 
coustic  temperature.  Then  all  nonacoustic  variations  of  temperature  in  I/, 
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are  part  of  the  turbulence.  For  example,  suppose  we  used  thermometers  to 
"measure"  a  "mean  temperature,"  and  found  these  readings  to  be  a  func¬ 
tion  of  height  only,  in  Vs;  typically,  they  would  decrease  with  height.  We 
could  model  this  in  terms  of  a  constant  ambient  T0,  the  average  over  height 
of  the  thermometer  readings,  plus  turbules  whose  amplitudes  AT  did  not 
cancel  on  average  and  depended  on  height  in  such  a  way  as  to  reproduce 
the  thermometric  results.  This  method  would  introduce  inhomogeneity  in 
the  turbulent  temperature  variations,  and  again  this  is  easy  to  do  with  a 
turbule  model.  These  ideas  may  be  worth  further  investigation.  However, 
this  approach  would  not  allow  for  the  curved  limiting  ray  paths  that  create 
shadow  zones  in  the  source  radiation  patterns.  This  curvature  is  small  and 
was  neglected  in  this  work,  because  we  treated  scattered  radiation  only. 

Another  approximation  that  was  made  in  the  theory  and  thus  in  the  de¬ 
veloped  code  was  that  the  source  and  detector  are  instantaneously  in  the 
radiation  zone  of  any  turbule,  i.e.,  that  RTS  and  RDT  are  large  compared  to 
the  larger  of  (A.  a2/ A),  at  all  times.  This  is  not  true  at  all  times  for  all  turbules 
in  realistic  scenarios,  in  which  a  few  of  the  turbules  may  actually  intercept 
the  source  and/or  detector  and  in  which  a  substantial  fraction  of  the  large 
turbules  spend  some  time  close  to  the  source  and  detector.  The  radiation 
zone  expressions  were  used  as  a  compromise,  because  they  yield  the  same 
results  as  the  conventional  theory  in  the  limit  of  (3  — >  0.  We  expect  that  cor¬ 
rections  to  this  approximation  would  not  alter  the  general  character  of  the 
total  |A(u>)|,  but  this  expectation  should  be  checked  in  the  future.  The  fully 
correct  theoretical  expressions  for  the  Born  approximation  scattered  waves 
in  near-field  regions  of  a  turbule  can  be  obtained  for  the  general  model 
used,  if  we  do  not  make  the  approximations  mentioned  just  before  equa¬ 
tions  (A-21),  (A-37),  and  (A-52)  and  just  after  equations  (A-38)  and  (A-39). 
But  then,  the  expressions  for  the  waves  are  extremely  complicated  and  dif¬ 
ficult  to  work  with,  so  much  so  that  it  seems  unlikely  a  three-dimensional 
model  could  be  handled. 

The  computer  code  Slow  Eddy  developed  in  this  work  is  not  optimized. 
It  can  model  a  three-dimensional  problem  only  roughly  on  a  PC  and  only 
because  the  new  time-shift  algorithm  is  so  much  faster  than  direct  compu¬ 
tation.  One  improvement  that  might  still  be  made  in  the  code  is  to  allow  the 
scattering  volume  to  be  different  for  each  turbule  size:  The  larger  the  tur¬ 
bules,  the  smaller  the  volume  needed.  But  a  recipe  for  choosing  the  optimal 
scattering  volumes  is  not  so  easy  to  attain,  because  the  choices  depend  on 
virtually  all  the  parameters  of  a  scenario.  It  is  probably  best  to  experiment 
with  this  code  to  help  determine  how  to  modify  it  to  use  scale-dependent 
scattering  volumes.  Unfortunately,  the  effective  volume  needed  increases 
inversely  with  turbule  size,  so  this  would  not  help  much  in  reducing  run¬ 
ning  time. 

Other  facets  of  an  experimental  scenario  are  that  source  waves  reflected 
from  the  ground  will  be  incident  on  the  turbulence  and  that  the  source 
and/or  the  detector  may  be  anisotropic.  The  latter  could  be  modeled  by 
suitable  modifications  of  the  computer  code,  along  with  some  further 
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theoretical  development.  The  former  would  require  substantially  more  the¬ 
oretical  development  and  code  alterations.  The  general  character  of  scat¬ 
tered  wave  spectra  may  be  altered  appreciably  by  both  effects. 

Overall,  the  frequency  spectra  predicted  by  this  model  of  time-dependent 
turbulence  seem  to  have  the  same  general  characteristics  as  some  observed 
spectra.  How  well  this  model  applies  to  many  different  sets  of  experimental 
results  needs  to  be  tested.  If  it  seems  to  work  fairly  well,  then  quite  likely,  it 
could  be  used  to  analyze  acoustical  remote  sensing  data  with  a  view  toward 
inferring  detailed  information  about  atmospheric  turbulence  and/or  wind 
speeds. 
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Appendix  A.  Theoretical  Derivations 


A-l  Geometry  and  Model 

See  figure  1  of  the  main  report.  In  the  derivation  that  follows — 
t  =  clock  time 

x.  y  =  horizontal  coordinates,  2  =  vertical  coordinate 
r  =  position  vector 

w  =  ambient  wind  velocity,  uniform,  constant 

T0  =  ambient  temperature 

p0  —  ambient  pressure 

pQ  =  ambient  mass  density 

Rs  =  s  =  source  location 

R°  =  d  =  detector  location 

R T(t)  =  RT  +  wt  =  b  +  wt  =  location  of  turbule  "center" 

RTS  =  R T(t)  -  Rs 
Rdt  '  =  Rd  —  R T(f) 

The  DS  or  ground-fixed  frame  is  the  "unbarred  "  frame  (D,  S  —  detector, 
source).  The  T  (barred)  frame  is  comoving  with  the  wind  ( T  =  turbule). 
Their  Cartesian  coordinates  are  related  by 


Xi  =  xt  —  Wit,  t  —  t  Galilean  transformation.  (A-l) 

The  total  (temperature,  pressure,  mass  density,  and  fluid  velocity)  fields  in 
the  DS  system  are 

T(r,  t)  =  T0  +  Tt(  r,  t)  +  T«(r,  t)  (A-2) 

p(r,  t)  =  p0  +  pt(r,t)  +  pa (r ,  t)  (A-3) 

p( r,  t)  =  p0  +  pt{ r,  t )  +  pa{r,  t)  (A-4) 

v(r,f)  =  w  +  vt(r,  t)  +  v0(r,  t)  ,  (A-5) 


where  the  subscripts  ( t .  a)  =  (turbulent,  acoustic).  In  the  T  system,  these 
fields  are 


T(v,i) 

—  T0 +  Tt{r,t) +  Ta{r,t) 

(A-6) 

=  Po  +Pt{r,t)  +pa{f,t) 

(A -7) 

v(? ,i) 

=  Vt(r,f)  +  V0(r,t) 

(A-8) 

Pi  r,t) 

=  Po+Pt(r,t)  +  pa(f.t)  , 

(A-9) 
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where 

T(r.F)  =  Ti(r.t) . v,(r./)  =  v,(r.  f)|r=r-w/  (A-10) 

and  similarly  for  acoustic  quantities. 

We  assume  that  the  source  emits  uniformly  in  4tt  and  that  the  ground  is 
totally  absorbing;  i.e.,  we  ignore  ground  reflections. 

A-2  Acoustic  Equations 

In  the  barred  (comoving)  frame,  the  relevant  acoustic  equations  are  the 
usual  ones: 

dtV„i  +  P~ld<pa  =  -  {Tt/p„T(,)  diVa  ~  Vu/i/'T,  ~  '  (A'11) 

dtp,,  +  ^hpJ),Vai  -  -Vjjdip,,  +  S„  .  (A-12) 

The  partial  derivative  symbols  are  dt  —  d/()t  and  0,  —  d/dx,  and  x,  — 
x.  x-2  =  y.  and  x:>,  =  c  in  the  stationary  frame.  The  F)r,  etc,  are  partial  deriva¬ 
tives  in  the  comoving  frame,  where  is  the  ratio  of  specific  heat  at  constant 
pressure  to  the  specific  heat  at  constant  volume. 

In  these  equations,  we  have  dropped  terms  in  pT,  since  this  is  second  order 
in  vT,  and  used  the  equation  of  state  to  put  pT/pa  +  Tt/T0  =  pT/p0  ~  0. 
Also,  we  have  included  a  scalar  acoustic  source  function  Sfl  to  be  chosen. 
Note  that  the  right-hand  sides  of  these  equations  involve  the  turbulent  tem¬ 
perature  and  velocity  fluctuations  TT(f.  t).  vTj(r.  /);  in  this  work,  we  regard 
them  as  stationary,  i.  e.,  Tt  (r),  vTi  (r)  only,  with  no  explicit  time  dependence 
in  the  barred  frame.  They  may  be  stochastic  functions.  Thus,  in  every  de¬ 
tail,  the  description  of  the  turbulence  and  acoustics  in  the  barred  frame  is 
subjectively  the  same  as  in  the  standard  description1  in  any  reference  frame 
with  no  wind. 

The  acoustic  quantities  in  equation  (A-ll)  and  (A-12),  the  sum  of  the  waves 
created  by  the  source  and  the  waves  scattered  by  the  turbulence,  are  written 
as 


Pa  —  Pac  "T  Pas  •  V ,,,  —  V„r j  +  l\,si  • 

in  the  Born  approximation  where  the  subscript  e  stands  for  "external,"  i.e., 
for  the  source  waves  that  act  as  externally  applied  fields  on  the  turbulence, 
and  the  subscript  s  stands  for  "scattered,"  representing  waves  scattered  by 
the  turbulence. 


A-3  Source  Wave 


The  "external"  acoustic  pressure  and  velocity  waves  in  the  barred  (T) 
frame,  generated  by  a  scalar  source,  satisfy  equations  (A-ll)  and  (A-12) 
with  (vTj.  Tt)  =0: 

dtp,,,  +  ipp/hvacj  =  S„(  r.  t.)  (A-13a) 

’V.  I.  Tatarskii,  The  Effects  of  the  Turbulent  Atmosphere  on  Wove  Propngntion,  Keter, 
Jerusalem  (1971). 
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(A-13b) 


dtVaeA  +  p0  diPae  —  0  • 

where  Sa( r.  f)  is  the  source  term.  From  equation  (A-13b),  we  may  put 

Vae  =  -»  pat  =  ~P<A$  -  (A-14) 

where  4»  is  a  "velocity  potential."  Then  from  equation  (A-13a), 

-p0d?$  +  %p0V2$  =  Sa  ■ 

Dividing  by  7 hp0  yields 

(V2  -  £$)$  =  -±zSa  •  C2  =  7,,Po/Po  •  (A-15) 

where  c  is  the  adiabatic  sound  speed. 

In  the  unbarred  frame,  the  monochromatic  isotropic  point  source  is  located 
at  r  =  Rs  =  s  and  has  frequency  vs  =  (jJs/2tt: 

Sa{ r,  f)  =  A'S(r  -  s)  exp(— iu>st)  +  c.c. ,  (A-16) 

where  A!  is  an  amplitude  and  <)(r)  is  the  Dirac  delta.  So,  using  equation 
(A-l), 

Sa(r ,  t)  =  A'S{r  -  s  +  wf)  exp(— iwsf)  +  c.c.  (A-l 7) 


The  Green  function  solution  of  equation  (A-15)  is 


A'  exp(— iu>st) 
P0C 2  47T 


—  s  +  wf  —  (3R)  exp(i/csf?) 


+  c.c. . 
(A-l  8) 


where  /3  =  w  /  c,  and 


R  =  r  -  r'. 


(A-19) 


We  note  that,  with  rj  =  r'  —  s  +  wf  —  /3R, 


8  ...  —  85  dr/j 

—d(f  -s  +  vt-PR)  =  ~  = 

B,5(r'  -s  +  wt-  PR)  =  = 


,  85  ' 

-c  (  —  pj 


(A-20) 


dsi 


Also,  we  neglect  (l  ./JR)  =  -TlJR2  compared  to  {l/R)  ikgdjl,  =  iksRi/R, 
since  we  will  have  kjR  >  1  always.  So  we  have  from  equations  (A-14),  (A- 
18),  and  (A-20) 


poe(r,  f)  =  p0c(iks  +  /3  •  Vs)  $  +  c.c. ,  (A-21) 


vQ«(r,  f)  =  (iks  +  /3  ■  Vs)  Ih  -1-  c.c. ,  (A-22) 
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where 


A1 

f(fr'^6( r'  -  s  +  wf  —  (3  R)  exp (ik„R)  exp(— iwst)  . 

(A-23) 

4?r f>„c-  j 

A' 

fd’r'Wdir'  —  s  +  wt  —  /3  R)  cxp(iksR)  exp . 

(A-24) 

4-/V2  J 

The  integrals  in  equations  (A-23)  and  (A-24)  can  be  done.  Transform  the 
variables  by 


Vi  =  r'i  ~  +  irjt  -  ;1jR  . 


(A-25) 


Then 


dm  -  ,  -9-  .  dft  r)R  ,  A  _  x 

+  Slnce  WJ  =  ~7F,  =  ~R>'  (A‘26) 


So  the  Jacobian  is 


,/(£)  =  i  +  0.r. 


Then  the  integrals  involve  S(tj),  so  we  get 


A!  exp(-iujst) 
4  7 rp0c2  R  +  /3-R 


exp(iksR). 


Ui  =  R&  . 


(A-27) 


(A-28) 


Here,  from  equation  (A-25) 

R  =  r  —  r'  =  r  -  [77  +  s  —  wf  +  /3R]v=q  =  R  —  (HR  ,  (A-29) 

where 

R  =  r  —  s  +  wf  =  r  —  s  (/)  (A-30) 

is  the  instantaneous  vector  from  the  source  to  the  observation  point  r.  Then 
taking  R  •  R  from  equation  (A-29)  yields 

Tf  =  (W)2  -  2/3  •  W~R  +  /327f\  or 
(1  -  ,32)  R2  +  2/3  •  Rs7?  -  (Wf  =  Q  . 

The  solution  is 

7?  -  V2  j-/3  •  R  +  ^(l  -  ;32)  (Wf  +  (/3  •  Ra)2| 

where 

'•f  =  (l  —  ;32)  1  .  as  in  relativity. 


(A-31) 

(A-32) 


Then 
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Thus 


OPar 

On 


A'k ;  RR, 

4~n‘2  (Dsf 


exp (—iajftt)  exp (ik^R)  +  c.c. 


(A-41) 


As  i  — •  0.  D’s  — *■  if .  R.  — »  R  \  —  1,  and  from  equation  (A-38),  we  get 

A'  ik ■  _ . 

Pn,  (  f  k)  =  — — =k  exp(-/^st)cxp(ik<I?s)  .  (A-42) 

47rr  yy 

Thus,  we  could  put  in  general 


-ik  A' 

— —  =  TV .  (A-43) 

C 

where  o'  =  length,  and  Vc  is  thus  the  pressure  amplitude  of  the  wave  from 
the  source  at  a  distance  R  —  o'  from  the  source  in  a  still  atmosphere.  This 
relation  identifies  the  significance  of  the  arbitrary  amplitude  A' . 


A-4  Scattered  Waves 

In  Born  approximation,  the  scattered  waves  satisfy  equations  (A-ll)  and 
(A-12)  with  ( pas .  va/ij)  on  the  left-hand  sides,  but  with  (pac,  vaei)  on  the  right- 
hand  sides.  By  combining  these  equations  again  to  eliminate  <9,  A/si  and  ap¬ 
plying  the  usual  solenoidal  condition  div n  =  0  on  the  turbulent  velocity, 
we  get  the  Monin  wave  equation 

(Of),  -  c2df)  pas  =  -T~ldi  (TjOrPar)  -  2pl)0,F)j  (vTjv„n)  ,  (A-44) 

where  two  terms  have  canceled  via  equation  (A-13a),  because  we  assume 
that  SQ  =  0.  where  the  turbulence  is  nonzero. 

We  write 


Pas  =  pis  +  Pas  ■  (A-45) 

where  pls  is  due  to  the  turbulent  temperature  fluctuations  Tr,  the  first  term 
on  the  right-hand  side  of  equation  (A-44),  and  p'af.  is  due  to  the  turbulent 
velocity  fluctuations  vTj,  the  second  term  on  the  right-hand  side.  Then  the 
Green's  function  solutions  of  equation  (A-44)  are 


Opac(f.t’) 


0*  i  ^  V =f—\r—r’\/c 

->  (r'.f-  If-r'I/c)  • 


(A-46) 

(A-47) 


where  r  is  now  the  observation  point.  After  doing  the  indicated  derivatives, 
we  will  put  r  =  d  -  w i.  d  =  R/;.  the  detector  location. 
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A-5  Temperature  Scattering 

From  equation  (A-41),  we  have 


r'.t') 


dr' 


A'kj  RRj 

i'=t-\f-r'\/c  47TC72  ( D  )3 


exp {iksR)  exy>(-iu!gt)  exp (iksR  )  +  c.c. 


(A-48) 


where 


but  now, 


TT”  -  -/ 

R  =  r  —  r 

Ri  =  R-  -  ;3,R .  i?  =  72  {  -/3  ■  Rs  +  2f } 


W 


=  RSJ(1-02)+(i3-R)2 


R*  =  r[  -  st  +  Wit,  -  0iR“  . 


-^d 


(A-49) 

(A-50) 

(A-51) 


(A-52) 


In  the  expression  for  p„s  of  equation  (A-46),  when  we  take  dj dru  we  again 
neglect  (d/dn)  (1  /Rd)  compared  to  ( iks/Rd)(d/dfi )  (R  +  Rd),  since  we  will 


have  ksRd  >>>  1  also. 


Similarly,  we  neglect  ( d/d)[ri  (RR,)  / (D)3]  also.  So  we  will  need 


_8_ 

dfi 

But 


exp 


ike 


R  +  Rd^j  j  =  iks-^r  [R  +  exp  [tfcs  [r  +  Rd)]  .  (A-53) 


8Rd__^d  dR_  3R_dR^  _  OR  j-d 
dfi  1  '  dfi  dRj  dfi  3  dR*  1 


(A-54) 


Using  equations  (A-50),  (A-51),  and  (A-54),  we  get 


dR 


dRs  D 


D° 


whereby 


d  -  -A  —d  / S  —  O  -\  R«  R 

—  (R  +  Rd)  =  Ri(D  -0-  R*  +  p2R)  -jjr  ■ 


Then 

d_ 

dfi 


(A-55) 


exp 


iks  (R  +  Rrf)]  =  iks  (r-R/D3^  exp  [iks  (R  +  Rd)]  .  (A-56) 


So  from  equations  (A-46),  (A-56),  (A-41),  and  (A-43),  we  get 
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(A-57) 


pL(r- *)  =  ~4~f~A::2  exi,I,/''s(/?  +  exp(-^0  +  c.c. 


Now,  the  turbule  is  localized  at  r'  =  b;  that  is, 


I'r(r,)  =  ^(f,-b)  =  7’(0 


(A-58) 


Define 


£  =  r  -  b  .  (A-59) 

and  change  variables  in  the  integrand.  Discard  the  £  terms  in  the  R.R,. 
etc,  in  the  numerator  and  denominator  of  equation  (A-57),  because  R(l  3> 
£.  R  >  £.  But  we  must  keep  terms  of  order  £  in  the  exponent.  That  is,  we 
write 


exp  iks  (^R  +  Rrl^j  rs  exp  iks  ^ 7?  +  Rrl^j  ^  exp(— iK  ■  £)  .  (A-60) 

and  this  defines  the  vector  K.  Then  we  put 


Then 


where 


R 


w t  =  d(f).  r;  =  b  +  £  . 

(A-61) 

-  b  -  w t-£  =  Rdt  -  £  . 

(A-62) 

Rd  =  Rdt  -  RDT  •  £  . 

(A-63) 

DI  =  d  —  b  —  w/  . 

(A-64) 

Now  we  need  R  as  a  Taylor  series  in  £.  Define  R()  by 

iHt)  =  a0  +  t-vbii. 


Here,  Rv  is  (i?)^_0  =  (R)f,=^,  with  R  defined  by  equations  (A-49)  to  (A-51). 
After  some  algebra,  the  previous  equations  yield 


<A-65> 

where  here 

^  =  ±-(su  +  R?t*j)K,-  (A-66) 
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with 


Rn  = 


RJS  -  i3jR 


DT 


y2{-/3-R  S  +  Ds} 

Ds  =  [(l  -  d2)  (i?s)2-  (/3  •  Rs) 
i?n,.  =  Rf  -  8,Rn  . 


(A-67) 

(A-68) 

(A-69) 

(A-70) 


Then  for  pjs  (r  =  d  -  w t..  t)  (at  the  detector),  we  get  from  equation  (A-57) 
^(detector)  =  AT(t.)  exp (-iu>st)  +  c.c. ,  (A-71) 

where  AT(t)  is  a  complex-valued  time-dependent  amplitude. 

k*  Tea'  R'&ADT 


A1  (t)  =  - 


exp [iks(R0  +  RDT)]T(K)  .  (A-72) 


47rT0  72  (D  )3Rdt 


with  K  given  by  equation  (A-65)  and  the  other  vectors  by  equations  (A-66) 
to  (A-70).  Here,  T( K)  is  the  Fourier  transform  of  T(^),  given  by 


T(K)  =  Jd3ZeM-iK-0T(1i).  (A-73) 

This  is  accurate  only  if  the  source  and  detector  are  in  the  far-field  region  of 
the  turbule  at  any  given  instant.  That  is,  we  must  have  RTS  >  A«,  Rts  >  a, 
and  Rts  3>  a2/\s,  where  a  is  the  scale  length  of  the  turbule,  As  =  2ir/ks, 
and  similarly  for  RDT. 

This  is  also  valid  to  all  orders  of  /3;  computationally,  we  should  use  it,  al¬ 
though  probably  only  first  order  in  j3  is  needed,  since  we  expect  8  <0.05 
in  most  realistic  situations,  except  in  strong  storms.  Note  that  (5  =  0.01 
yields  w  =  3.4  m/s  ss  7.6  mi/hr,  so  f3  <  0.05  — »  w  <  38  mi/hr,  reason¬ 
able.  But  terms  of  order  01  in  K  and  in  exp[iks(R(t  +  RDT)]  might  shift 
the  phases  noticeably,  since  ks  =  us/c  =  (27r)(500)/340  =  9.24  m_1.  It 
would  be  accurate  enough  to  keep  only  first-order  terms  in  (3  in  the  term 
,R0(R0  •  R dt)/(Ds)3Rdt,  but  that  would  save  very  little  time  in  numerical 
computation,  so  the  higher  order  terms  are  kept. 

For  simplicity,  we  use  a  Gaussian  turbule, 

T(0  =  (AT)exp(-e2/«2)-  (A-74) 

The  size  or  scale  of  this  is  a.  Then  the  Fourier  transform  is 

f  (K)  =  7r3/2a3(AT)  exp(-A2a2/4)  .  (A-75) 
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effective  incident  fj  (0)  spherical  outgoing 

plane  wave  amplitude  scattering  amplitude  wave  from  scatterer 


where  we  have  defined  RDT  =  r.  RDT  =  r.  R/  s  =  ks.  and,  for  ft  —  0. 

7?0  =  kRTs.  D'  =  Rrs.  R(1  •  Rpr  =  cos  0  . 

with  (9  the  scattering  angle.  This  A1  (t)  agrees  with  equation  (2)  found  in 
an  earlier  paper  by  Goedecke  and  Auvermann;2  it  is  the  Monin  result  for 
acoustic  scattering  by  stationary  turbulent  temperature  inhomogeneities. 

The  AT  is  the  temperature  fluctuation  amplitude.  For  isotropic  homoge¬ 
neous  turbulence,  it  scales  with  a  as  |AT„|  =  j AT,  |  (a /a , ) 1  \  as  discussed 
in  the  main  text.  If  Tu  is  the  mean  temperature  T,  then  half  the  A7),  must 
be  taken  negative  and  the  other  half  positive,  but  with  the  same  magnitude 
for  each  a. 


A-6  Velocity  Scattering 

From  equations  (A-39),  (A-52),  and  (A-47),  we  see  that  we  need 

+*,)). 

Previously,  we  found  in  equation  (A-56) 

A  exp [iMB  +  =  ik,  (r1,  R/D ^  exp [iks(R  +  /?,,)]  . 

We  obtain  the  derivative  of  this  here  by  once  again  neglecting  OR/ Of ,  com¬ 
pared  to  k(<R  OR /Of,,  valid  since  ksR  »>  1.  So  we  get 
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Now  we  put 

%(f')  =  Vj  (r/  -  b)  =  Vj  (£)  ,  (A-79) 

and  we  put  £  =  0  in  the  terms  in  the  square  bracket  in  the  integrand.  So 

there  R,  =  7?f  .  7?,  =  R0i .  i?  =  i?.0 .  etc,  as  in  equations  (A-66)  to  (A-70), 
and  also 

exp[i/,;s(f?.  +  Rd)}  =  exp[iks{Ro  +  RDT )]  exp(-iK  •  £)  (A-80) 

as  in  equation  (A-60),  where  K  is  given  by  equation  (A-65).  Therefore,  for 
r  =  d  —  w  f, 


Pas  (detector) 


kiv,a'  (R„f  (R,ROT)  (fiDT  -m) 


|  27 rc  72(Z?sj4  RDT 

x  exp[i/cs(^o  +  jRDT)]  |  exp(— iujsi)  +  c.c. 

=  Av(t)  exp[—  iu)st]  +  c.c. . 

where  Av(t)  is  defined  by  the  brace.  Here, 

*(K)  =  J  d3^expHK-0v(0  • 

For  simplicity,  we  use  a  Gaussian  turbule 

v  (£)  =  (fi  x  ()  exp(— £2/ a2) 
as  in  Goedecke  and  Auvermann's  paper.  Thus 


f(K)  =  ftxiVx  y*d3^exp(— iK  ■  £)  exp(-£2/a2) 


ia37r3^2H  x  Vk  exp(— iT2a2/4)  or 


v(K)  = 


ia57r3/20 


(A-81) 


(A-82) 


(A-83) 


(n  x  k)  exp(— A'2a2/4)  .  (A-84) 


For  isotropic  homogeneous  turbulence,  the  energy  transfer  rate  per  unit 
mass,  e,  from  one  turbule  size  to  the  next  satisfies 

e  c x  (z>2)  ( v/a )  =  (independent  of  a)  fa  ft3a2  . 

Thus  Q(a)  scales  with  a  as 

/  \  2/3 

0(a)  =  Q(aJ  f  —  j  ,  aj  =  amax  =  largest  turbule  scale  length  used.  (A-85) 
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So  we  put 


•  5_3/2 

*(K)  =  )«7/3)  (ft  x  k)  exp(-A'V/4)  •  (A-86) 

In  equation  (A-84),  ft  is  a  "randomly"  oriented  unit  vector.  If  the  turbulence 
is  assumed  to  be  isotropic,  then  ft  has  a  uniform  probability  distribution  in 
47 r;  that  is,  if  we  write 


ft 


=  e;  sin  Bv  cos  ov  +  e_,  sin  Bv  sin  6U  +  e;i  cos  0V  . 


then 


PA0JM» 

P0  (A>  )  (l0(! 


sill  0i>(10{} 


C,(«u) 


P.(A,)  = 


sin 

2 


where  (Pu  (0it) .  P  (os>))  are  the  probability  distributions  of  (0il.(f>il).  If  the 
turbulence  is  not  taken  to  be  isotropic,  then  a  spherically  symmetric  ro¬ 
tating  turbule  model  may  not  be  adequate,  and  also  ft  is  not  uniformly 
distributed  in  An. 


Consider  the  result  equation  (A-81)  in  the  limit  /3  =  0.  Insert  b  =  0  =  b. 
Then  RT5  =  k:  RDT  =  r  =  d:  R0  =  RTS;  TP  =  RTS:  RDT  =  r:  so  we  get 


Vp(l  exp (iksRTS) 


effective  incident 
plane  wave  amplitude 


RTS 


irz  cos0(r- v(K)) 

b - v - " 

fv  ((*■<!>) 

scattering  amplitude 


exp  (rksT) 
r 


(A-87) 


spherical  outgoing 
wave  from  scatterer 


where  here 

K  =  (f  -  k)  • 


(A-88) 


Note  that  K  •  v(K)  =  0,  since  Vc  •  v  (^)  =  0.  whereby  r  •  v(K)  =  k  •  v(K) 
for  (3  =  0.  This  scattering  amplitude  is  the  same  as  equation  (2)  of  Goedecke 
and  Auvermann's  paper.  (This  is  the  Monin  result  for  scattering  by  turbu¬ 
lent  velocity  fluctuations.) 

We  note  that  K  •  v(K)  =  0  in  general,  for  any  (3.  for  the  rotating  sphere 
model  (see  eq  (A-82)).  In  fact,  this  must  be  true  for  any  turbulent  eddy  that 
satisfies  V,  •  v  (£)  =  0.  because  then 

j<P£  exp(— v'K  •  OV,  •  v  (0  =  0  =  /K  •  ?(K)  .  (A-89) 

where  integration  by  parts  and  discard  of  the  surface  integral  at  "spatial 
infinity"  yield  the  second  equality.  This  surface  integral  must  vanish  since 
each  turbule  is  localized.  (We  assume  that  turbules  are  "wholly  contained" 
in  the  scattering  volume.) 
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Appendix  B.  Computer  Algorithm 


B-l  Detector  and  Source  Locations  and  Wind 

The  calculations  are  done  in  the  "primed"  coordinate  system,  which  is 
ground-fixed  but  in  which  the  wind  is  along  the  x'  axis, 

WXI  —  w  >  0.  Wy>  =  W-l  =  0, 


as  in  figure  B-l.  Define 


a:CM  =  +  Xs). 

vCM  =  \(vD  +  ys)- 

The  coordinate  transformation  is  then 

J  _  _ 

x'  =  (x  —  a:CM)  cos  a  +  (y  —  yCM)  sin  a,  (B-l) 

y1  =  —  (x  —  a:CM)  sin  a  +  (y  —  yCM)  cos  a  , 

and  components  of  all  position  vectors  are  related  in  the  same  way.  The 
user  inputs  the  Cartesian  components  of  (R°,  Rs)  and  the  magnitude  w 
and  azimuthal  angle  a  of  the  horizontal  wind  velocity  w,  in  the  "original" 
(unprimed)  frame.  The  algorithm  calculates  the  primed  components  of  all 
vectors. 


Figure  B-l.  "Primed 
coordinate  system. 
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B-2  Spectral  Analysis 

Consider  the  Fourier  transform  of  equation  (8)  from  the  main  report, 

roc 

Pa  si'*'')  —  /  r//exp(/o,'/)  [.47  (/)exp(— D,'.,f)  +  AT*(/)  oxp(/a,\,/)] 

•/  —  !X 


or 


Pas ~  A^  (a.’  —  a.’*)  +  A1  ’  (a,’  +  Uls)  ■ 

Therefore,  all  we  need  to  calculate  is 

A7  (a;)  =  f  (It  exp(iujt)AT  (t) 

./  —  DC 


(B-2) 


(B-3) 


and  then  use  equation  (B-2)  if  we  wish.  But  the  spectral  broadening  and 
shifting  is  all  contained  in  .4 7  (v),  and  similarly  for  the  relation  of  f>„.s(u/)  to 

A1 


Since  we  use  fast  Fourier  transform  algorithms,  we  must  establish  the  min¬ 
imum  number  of  time-series  points  needed  for  adequate  accuracy.  First, 
equation  (42)  from  the  main  report  shows  that  the  frequencies  in  ATW) 


are  given  (approximately)  by  vs(3  ■  (kD7  (t)  -  R TS(t)'j  ■  This  varies  from 
a  minimum  of  -2 f)us  to  a  maximum  of  2fli /H;  so  from  the  Nyquist  criterion, 
we  need  for  AT (ui) 


/-'max  =  ((2 A t)~l)  >  2ftus  -  At  <  (4 ftus)-1  .  (B-4) 


For  example,  for  us  =  500  Hz,  0  =  0.01,  we  need  At  <  0.05  s  —  //max  > 
10  Hz.  The  computer  code  uses  an  algorithm  that  yields  these  numbers  for 
=  500  Hz,  0  =  0.01.  That  is,  it  chooses 


"max  =  integer  value  >  2 ftvs  .  At  =  (2//max)  1  .  (B-5) 


Second,  we  need  appropriate  frequency  resolution.  We  can  obtain  an  esti¬ 
mate  of  Tr  (the  minimum  total  time  duration  of  a  measurement  that  en¬ 
sures  adequate  frequency  resolution)  by  examining  the  scattered  wave  am¬ 
plitude  associated  with  large  turbules  with  ksa  >  1.  Consider  figure  B-2. 
For  what  values  .Tj  of  .r  is  the  Gaussian  reduced  to  perhaps  c- 1  of  its  maxi¬ 
mum  value?  If  kfO  1.  .r  will  be  close  to  D/2.  A  little  algebra  yields 


D  1  (D-  +  4  Irf!'2 

2  y/2k„a  Dli 


(B-6) 


Now  use  equation  (42)  from  the  main  report  to  obtain  //,  when  .r  =  D/2 
and  :r  =  .  We  easily  get 


",  =  Pr=D/2  ~ 


\pl  ir 

770 


f  D~  +  Ah'1  \ 

V  2  Dh  ) 


(B-7) 
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Figure  B-2.  Estimating 
minimum  measurement 
time  duration. 


Essentially,  only  —v1  <  v  <  will  be  in  the  Fourier  transforms  of  the 
amplitudes  A(t). 

As  w/a  gets  smaller,  so  does  u1.  The  quantity  D 2  +  4/?2/2D/i  has  minimum 
value  2  for  h  =  D/2 .  So 


)min 


(B-8) 


Adequate  frequency  resolution  requires  perhaps  Au  =  TR  1  <  10  1  (i^  )min. 
We  suppose  that  amax  —  10  m,  u)min  =  1  m/ s.  Then 


m  \  1  r>  /  /  \  107ramax 

>  10/  (^)min  =  — - »  200  s  . 


y/2- 


w„ 


(B-9) 


Independently  of  this,  we  would  like  the  total  time  T  to  be  long  enough  to 
contain  almost  all  the  scattered  pulse  from  one  turbule  that  passes  "sym¬ 
metrically"  over  the  source /detector.  For  small  turbules,  the  scattered  pulse 
amplitudes  are  approximately  proportional  to  (RDTRTS)  1 .  For  source  and 
detector  positions  both  near  the  ground  and  a  turbule  passing  directly  over¬ 
head,  as  in  figure  B-3,  the  scattered  amplitude  is  essentially  a  maximum 
when  the  turbule  is  at  position  2,  proportional  to  where 

Rmin=((D/2)2  +  h2f2  .  (B-10) 


At  positions  1  and  3,  the  amplitude  is  proportional  to  R  2,  where 

R2  =  {(R1+  D/2)2  +  h 2)  5  ((^  +  D)2  +  h 2)  *  .  (B-ll) 

For  most  scenarios,  for  turbules  of  interest,  so  we  could  neglect  h  in 

equations  (B-10)  and  (B-ll).  Then  we  could  require  that  we  limit  the  length 
of  the  scattering  volume  along  the  wind  direction  (which  is  along  +x  in  this 
simple  analysis)  by  requiring  R2  >  16f?1^lin;  this  would  include  all  signals 
within  6  percent  of  the  maximum  from  a  turbule  that  starts  at  x  =  —  R1  and 
travels  the  total  distance  2R^  +  D  in  time  T.  For  h  «  0,  this  yields 

(Rt  +D)>2D^R,>D.  (B-12) 
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Figure  B-3.  Estimating 
maximum  measure¬ 
ment  time  duration. 


Then,  the  total  length  L  of  a  tube  along  the  wind  would  satisfy 

L  =  D  +  2/?,  >3D  =  3|R/>S  |  .  (B-13) 

The  computer  algorithm  uses  the  equality,  and  then  chooses 

T  =  L/ir  =  3|R/w|/«- .  (B-14) 


or  uses  T  =  Tu  given  by  equation  (B-9),  whichever  is  larger.  For  most  cases 
of  interest,  equation  (B-14)  yields  the  larger  value. 

The  minimum  number  of  points  N  needed  in  the  discrete  time  series  for 
the  A(t)  is  then  obtained  by  the  use  of  equations  (B-4)  and  (B-14), 


N  >  T/At 


3Rns4,3vs  6 ,  _ 
—  A  c 


ns 


w 


(B-15) 


or  by  the  use  of  equations  (B-4)  and  (B-9), 


Ar^^/A,  107romax2pns  10 
A  s  r/A' = 


(B-16) 


whichever  is  larger.  Note  that  usually  BDS  »  o.max. 

Since  A'  is  used  to  dimension  arrays,  the  FORTRAN  source  code  actually 
defines  a  default  value  N  =  6048  as  a  parameter.  The  user  should  verify 
that  this  is  large  enough  to  satisfy  equations  (B-15)  and  (B-16),  and  if  it  is 
not,  change  it  and  recompile.  The  values  of  T  and  L  are  calculated  inter¬ 
nally  by  putting  (T  =  (N  —  1)A/.  L  =  wT)  . 


B-3  Turbule  Model  Parameters 

As  shown  by  Goedecke  and  Auvermann,1  equations  (67)  and  (68),  for  Kol¬ 
mogorov  spectra  of  isotropic  homogeneous  turbulence,  the  structure  pa¬ 
rameters  are  related  to  the  turbule  parameters  by 

C~.  =  (3.78///)  (4>)  (A  TJa^f  J8/3  .  (B-17) 

C;2  =  (0.68///.)  (</>)  (jljf/2/'*)2  Jj.,/3  . 

]G.  H.  Goedecke  and  H.  J.  Auvermann,  "Acoustic  scattering  by  atmospheric  turbules," 
].  Aconst .  Sac.  Am.  102  (1997),  pp  759-771 . 
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where  successively  smaller  turbule  sizes  aa  are  related  by 


",,+lK  =  e  =  NrJ  1 


(B-18) 


where  a ,  =  «max  is  the  outer  scale  length,  <f>  =  n,a,  =  naaa  is  the  packing 
fraction  (the  same  for  all  sizes),  and  nQ  is  the  number  density  of  turbules  of 
size  a.  ( AT, .  0. , )  are  the  (temperature  fluctuations,  angular  velocity)  of  the 
largest  turbules,  and  here 


J,  =  /  dx  xs e  T  =  -T 


lx,  /«+  1 


Jo 


(B-19) 


where  T  is  the  gamma  function.  From  tables  of  F  functions. 


l  , 

/II 

II 

co 

V6 

1^  i 

( 17 

hi/z  — 

Vi 

v  6 

0.47. 


=  — r  (— 
12  v  6 


(B-20) 


0.S 


The  (A Tq,Qq)  for  size  aa  are  scaled  from  (AT,,!},)  as  in  equations  (31) 
and  (32)  from  the  main  report.  The  parameter  /x  =  In  Nj  in  equations 
(B-14)  and  (B-15)  is  calculated  after  the  inner  scale  length  ajvs,  the  small¬ 
est  turbule  size,  and  the  number  of  sizes  Ns  in  the  cascade  are  chosen  by 
the  user: 


aivs/°-i  =  exp[-/j,(Ns  -  1)]  — >•  y  =  ~(NS  -  1)  1  ln(ajvs/aJ  .  (B-21) 


Default  values  are  Ns  =  6,  a1  =  1.5  m,  a^s  —  a6  =  4.6875  cm,  which 
yield  /i  =  In  2,  Nj  =  2.  i.  e.,  a  halving  of  sizes  with  each  step.  Doubling 
Ns  while  keeping  (ax ,  )  unchanged  should  approximately  double  the 

running  time  of  the  code.  The  effects  of  changing  Ns,  air  and/or  a^.,  had 
not  been  investigated  at  the  time  of  this  writing. 

Values  of  (C'*.C%)  may  be  input,  in  meter-kilogram-second  units,  and  then 
the  appropriate  values  of  (Tt  =  AT3  /T0 ,  rt,  =  £l1al/ c )  are  calculated.  If  the 
values  =  0  are  input,  then  the  code  sets  default  values  of  ( AT) .  ff  a ) 

by  putting 


AT,  /T0  =  rT,  fi1  a,  /c  =  rt, ,  (B-22) 

where  (rT,r„)  may  be  input.  The  algorithm  begins  calculating  with  the 
largest  turbules  (size  a,).  It  places  the  first  tube  of  turbules  of  square  cross 
section  d1  =  a,  / 01/3  along  the  wind  direction  (T-axis)  just  on  the  ?/  side 
of  the  origin  and  centered  along  the  x'-axis.  The  origin  is  chosen  halfway 
between  the  source  and  the  detector.  That  is,  the  centerline  of  the  tube  is 
y’  =  dj 2,  £  =  Rs( 3)  +  d0  +  d,/2.  (The  base  of  this  tube  is  calculated  to  be 
d0  m  above  the  height  of  the  source  or  detector  above  ground,  whichever  is 
greater.)  The  length  of  the  tube  is  chosen  by  equation  (B-13).  Then  at  t  =  0, 
one  turbule  is  placed  randomly  inside  each  cubical  cell  of  side  length  d1  in 
the  tube,  such  that  each  turbule  center  can  have  (x',  y\  z’)  displacements 
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from  its  cell  center  from  0  to  ±d/ 4  with  equal  independent  probabilities. 
Then  ( AT(t ).  Ar(t))  are  calculated  and  accumulated  for  all  these  turbules. 
Any  part  of  a  turbule  that  leaves  the  tube  (at  x'  =  L/2)  during  time  T  is 
introduced  at  x'  =  —L/ 2  +  irAt.  the  next  time  step,  thus  keeping  the 
frozen  turbulence  model  in  steady  state. 

Successive  tubes  for  the  largest  turbules  are  chosen  in  such  a  way  that  tubes 
chosen  later  are  further  from  the  origin,  in  general.  For  example,  the  sec¬ 
ond  tube  would  have  a  centerline  (//  =  — c/,  /2.  z'  =  Rs (3)  +  <7(l  +  <7,  /2).  The 
third  and  fourth  tubes  would  be  above  the  first  and  second,  respectively, 
with  centerlines  (;</  =  ±dj2.  z'  =  i?'s(3)  +  r/0  +  3<7,/2).  The  fifth  and 
sixth  tubes,  if  used,  would  be  one  step  further  away  in  ?/,  with  centerlines 
y'  =  ±MJ2.z'  =  Rs( 3)  +  r/0  +  dj2.  This  selection  continues  until  input 
limits  or  default  limits  are  reached. 

For  smaller  turbules,  successive  tubes  are  chosen  that  fill  all  the  largest 
tubes  used.  This  implies,  for  example,  that  if  a  ,  =  nj 2.  then  there  will 
be  four  times  as  many  tubes  of  size  d2  =  o2/d1//,i  =  dj 2  as  there  are  for 
size  d, .  It  is  impractical  to  calculate  turbule  by  turbule  for  the  smaller  tur¬ 
bules,  because  there  are  so  many  of  them.  Instead,  the  time-shift  algorithm 
is  employed  for  turbules  with  ( RTS .  RDT)  >  (A’.,r/2  or  k~l),  whichever  is 
larger,  as  described  in  the  main  report.  The  crucial  values  of  R1S  and  Rm 
for  this  determination  are  the  values  for  which  the  maximum  signal  occurs. 


Appendix  C.  Quick  Doppler  Formula 


Refer  to  figure  1  of  the  main  report  (p  3)  for  the  geometry.  The  turbule  is 
at  R T(t)  =  b(f)  =  b0  +  w t.  where  w  —  wind  velocity.  Pulses  are  emitted 
from  the  source  S  located  at  position  R5  =  s  at  times  tf ,  tf  >  if  with 
if  -  u>s  =  source  angular  frequency.  Then  the  distances  from  the 

source  to  the  turbule  at  times  (if.  if)  are 


RjS  =  |b(if)  -s|  =  |b0  -  s  + wif|  , 
RTS  =  |bo  -  s  +  wtf  |  =  \KTS  +  w^' 

=  R.fs  +  w  .  R f  »(tf  -  if) . 


w  (if -if) 


Thus  these  pulses  get  to  the  turbule  at  times  (if,  if)  given  by 


fT  =  tS 


1  c(  l  +  /3-Rf5) 


■dTS 

if  +  -M1-0  -Rf5) 

1  n 


r>TS  dT5 

+T  _  +S  ,  n2  _  +S  ,  ^2  m  *  •8T5^ 

*■  ~K  +  ^rW)~ 2  ~  2  ] 


if  +  \  Ks  +  w  •  - 1; )  J  (1  -  p  ■  r;»)  , 


TSuT  J.T 


K  ~K  =n~n  +  /9-R-X  (*;  -tO*(£-ff)(i  +  0-Rx  )•  (C-3) 

The  distances  from  the  turbule  to  a  detector  located  at  RD  =  d  at  times 

(if,  if) 


=  |d  -  b0  -  wt,  | , 

=  |d  -  b0  -  w/f  -  w(if  -  if)  I  «  Rf r  (w  •  Rf r)  (if  -  if )  ,  (C-4) 

so  the  times  (if,  if)  at  which  the  scattered  pulses  are  received  by  the  de¬ 
tector  D  are 


i  +  -Rf 4 
c 

if  +  KDT 


t!  +  -Rf 


d-Rfr  (if-if) 
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Thus  to  order  3, 

t?  -t?  =  (e  -?)( 1-0-  Rfr)  =  (t*  -  /f)[l  -  /3  •  (Rr  -  RD)  •  (C-6) 

Therefore  the  frequency  observed  at  the  detector  is 

«  ^(1  +  0  •  (Rf,r  -  R™)) .  (C-7) 

This  is  Wobs  at  D  at  time  =  tj  +  £|d  -  b„  —  w t[  |  ~  t*  +  j  (ft(s  +  R™  ) . 
Both  R,y  and  Rrs  should  be  evaluated  at  time  t[  in  this  formula  for  u\*k. 
To  first  order  in  / 3 ,  this  is  the  same  result  obtained  in  the  main  report  directly 
from  the  expressions  for  the  waves  (p]is(t).  pj,s(0)  at  the  detector. 
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Appendix  D.  Computer  Codes 


This  appendix  contains  descriptions  of  four  FORTRAN  codes.  They  are 
written  in  structured  FORTRAN  77  standard  form.  The  Microsoft  Inter¬ 
national  Mathematical  and  Statistical  Library  is  used  for  random  number 
generation  and  calculation  of  fast  Fourier  transforms.  The  Microsoft  portlib 
routine  TIMERQ  was  used  to  time  the  codes.  The  codes  are  input.for  used 
with  sloweddy.for  and  inputl.for  used  with  slowone.for. 

The  input.for  and  inputl.for  codes  write  parameters  to  files.  The  parame¬ 
ters  in  the  files  can  be  edited  by  a  text  editor  without  having  to  rerun  the 
input  codes. 

The  description  and  code  (with  a  list  of  symbols)  for  program  input  follow. 
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★ 

* 


* 

* 

* 


CHARACTER  CONFIG 
INTEGER  IT, NS, NT 

REAL  ALPHA, AP , A1 , C , D ( 3 ) , GAMT , GAMV , MW , NUS , PF 
REAL  S ( 3 ) ,YC1{20) , ZB , ZC1 ( 20 ) 

OPEN ( 1 , FILE= ' input . txt ' ) 

c  The  following  inputs  have  to  do  with  the  environment  for  the  turbules. 
PRINT  * 

PRINT  *,  'Enter  the  sound  speed  in  meters  per  second:' 

READ  *,  C 

PRINT  * 

PRINT  *,  'Enter  the  sound  frequency  of  the  source  in  Hertz:' 

READ  *,  NUS 

PRINT  * 

PRINT  *,  'Enter  the  reference  distance' 

PRINT  *,  'from  the  source  in  meters:' 

PRINT  *,  ' (If  not  sure,  enter  1)' 

READ  *,  AP 

PRINT  * 

PRINT  *,  'Enter  the  x,y,z  coordinates' 

PRINT  *,  'of  the  source  in  meters:' 

READ  * ,  S 

PRINT  * 

PRINT  *,  'Enter  the  x,y,z  coordinates' 

PRINT  *,  'of  the  detector  in  meters:' 

READ  * ,  D 

PRINT  * 

PRINT  *,  'Enter  the  boundary  layer  height  in  meters:' 

PRINT  *,  ' (higher  than  either  source  or  detector)' 

READ  * ,  ZB 

IF( (ZB.LE.S(3) ) .OR. (ZB . LE . D ( 3 ) ) )  THEN 

PRINT  *,  'The  turbule(s)  will  pass  through  source  or  detector.' 

PRINT  *,  'Please  run  the  input  program  again,  and  enter  a' 

PRINT  *,  'boundary  layer  height  greater  than  the  source  or' 

PRINT  *,  'detector  height.' 

STOP 
END  IF 


*  Program  Input 

*  Version  1.0:  11-21-98 

******************************************************************* 

*  This  program  enters  the  environment  parameters,  and  the  turbule 

*  parameters  needed  for  the  Slow  Eddy  program.  These  parameters 

*  are  written  to  the  file  "input . txt . " 
***************************************************************** 

PROGRAM  INPUT 
IMPLICIT  NONE 
INTEGER  NMAX 
PARAMETER ( NMAX- 6400) 
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PRINT  * 

PRINT  *,  'Enter  the  wind  speed  in  meters  per  second:' 

READ  *,  MW 

PRINT  * 

PRINT  *,  'Enter  the  wind  direction  in  degrees' 

PRINT  *,  'measured  counter  clockwise  from  the  x-axis:' 

READ  *,  ALPHA 

The  following  inputs  have  to  do  with  the  turbules  themselves. 
PRINT  * 

PRINT  *,  'Enter  the  size  of  the  largest  turbule  in  meters:' 
READ  *,  A1 

PRINT  * 

PRINT  *,  'Enter  the  packing  fraction:' 

READ  *,  PF 

PRINT  * 

PRINT  *,  'Enter  the  temperature  scale  parameter:' 

READ  *,  GAMT 


PRINT  * 

PRINT  *,  'Enter  the  velocity  scale  parameter:' 

READ  *,  GAMV 
PRINT  * 

PRINT  *,  'Diagrams  for  configurations  are  given  below:' 
PRINT  * ,  ' (hit  "enter"  to  see  more) ' 

READ  * 

PRINT  * 

PRINT  *,  '  Configuration  choices  ' 

PRINT  * 

PRINT  * 

PRINT  *,  '  - 

PRINT  *,  'a:  |  | 

PRINT  *,  '  - 

PRINT  * 

PRINT  * 

PRINT  -  - 

PRINT  *,  'b:  I  I  I  I 

PRINT  *,  '  -  - 

PRINT  * 

PRINT  *,  ' (hit  "enter"  to  see  more)' 

READ  * 

PRINT  * 


PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 


c  : 


PRINT  * 


PRINT  * 


PRINT  *,  '  -  - 

PRINT  *,  'd:  |  |  |  | 

PRINT  *,  '  -  - 


PRINT  * 


PRINT  *,  ' (hit  "enter"  to  see  more)' 


READ  * 

PRINT  * 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  *,  ' e: 

PRINT  *,  ' 

PRINT  * 

PRINT  * 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  * ,  ' 

PRINT  *,  'f: 

PRINT  *,  ' 

PRINT  * 

PRINT  * 

PRINT  *,  'g: 

PRINT  * 


User  defined  geometry:' 


PRINT  * 

PRINT  *,  'Enter  the  letter  (a-g)  corresponding  to  the' 
PRINT  *,  'the  appropriate  configuration  desired.' 

PRINT  *,  ' (make  sure  letter  is  lower  case)' 

READ  *,  CONFIG 

PRINT  *,  'Configuration  chosen  was  ', CONFIG 


PRINT  * 

PRINT  *,  'Enter  the  total  number  of  turbule  sizes:' 
READ  *,  NS 


c  If  configuration  g  is  chosen,  the  user  will  need  to  enter 
c  the  number  a  large  tubes,  and  the  y,z  coordinates  to  the 
c  center  lines  of  these  tubes. 

IF (CONFIG. EQ. 'g' )  THEN 
PRINT  * 

PRINT  *,  'Enter  the  number  of  large  tubes:' 

PRINT  *,  ' (maximum  number  allowed  is  20)' 

READ  *,  NT 

DO  5  IT=1 , NT 
PRINT  * 

PRINT  *,  'Enter  the  Y  and  Z  coordinates  of  the  center' 
PRINT  *,  'of  large  tube  number' , IT 
PRINT  *,  'in  dimensionless  units:' 

READ  *,  YC1 (IT) , ZC1 (IT) 

5  CONTINUE 

END  IF 


c  All  parameters  entered  above  are  written  to  the  file  'input.txt' 
c  used  in  the  Slow  Eddy  FORTRAN  program. 

WRITE (1,*)  C 
WRITE (1,*)  NUS 
WRITE (1,*)  AP 
WRITE (1,*)  S 

WRITE (1,*)  D 
WRITE (1,*)  ZB 

WRITE (1,*)  MW 
WRITE (1,*)  ALPHA 
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WRITE (1,*)  A1 
WRITE (1,*)  PF 
WRITE (1,*)  GAMT 
WRITE (1,*)  GAMV 
WRITE (1,*)  CONFIG 

WRITE (1,*)  NS 
WRITE (1,*)  NT 

DO  20  IT=1 , 20 

WRITE (1,*)  YC1 (IT) , ZC1 (IT) 
20  CONTINUE 

CLOSE ( 1 ) 


END 

* _ * 

*  -  List  of  Symbols  - * 

*  _ _ _ * 

* 

*  ALPHA  :  real;  angle  of  wind  relative  to  non-rotated  x-axis. 

*  AP  :  real;  'a  prime',  reference  distance  from  the  source. 

*  A1  :  real;  size  of  largest  turbule . 

*  C  :  real;  speed  of  sound. 

*  CONFIG  :  character;  abbreviation  for  'configuration'  used  to  choose 

*  predefined  tube  configurations. 

*  D(3)  :  real;  position  vector  to  the  detector. 

*  GAMT  :  real;  'capital  gamma  sub  T' . 

*  GAMV  :  real;  'capital  gamma  sub  v' . 

*  IT  :  integer;  index  for  largest  tube,  with  range  from  1  to  NT . 

*  MW  :  real;  'magnitude  of  wind' . 

*  NMAX  :  integer;  maximum  number  of  time  (and  frequency)  data  points 

*  NS  :  integer;  number  of  sizes  of  turbules. 

*  NT  :  integer;  number  of  largest  tubes. 

*  NUS  :  real;  'Greek  letter  nu  sub  s',  the  frequency 

*  of  the  source . 

*  PF  :  real;  'packing  fraction'. 

*  S(3)  :  real;  position  vector  to  the  source. 

*  YC1 (20)  :  real;  y  coordinate  to  center  line  of  largest  tube,  dimensionless. 

*  ZB  :  real;  height  above  ground  where  scattering  volume  begins. 

*  ZC1(20)  :  real;  z  coordinate  to  center  line  of  largest  tube,  dimensionless. 


Descriptions  of  Program  Slow  Eddy  and  its  eight  subroutines  and  a  list  of  symbols  follows 


******************************************************************************* 

*  Program  Slow  Eddy  * 

*  Version  1.0:  11-21-98  * 

******************************************************************************* 

*  This  program  calculates  the  pressure  amplitudes  as  a  function  of  time,  * 

*  and  their  Fourier  transforms  as  a  function  of  frequency.  This  program  * 

*  uses  input  from  the  file  'input.txt'  that  is  generated  by  running  the  * 

*  'input. for'  Fortran  program.  The  Slow  Eddy  program  is  designed  to  populate  * 

*  the  scattering  volume  with  tubes  of  turbules.  There  are  six  predefined  * 

*  configurations  (a-f)  and  one  configuration  choice  specified  by  the  user  (g).* 

*  The  tubes  are  populated  with  NS  different  sizes  of  turbules,  keeping  the  * 

*  packing  fraction  constant.  Each  smaller  size  is  obtained  by  halving  the  * 

*  previous  size.  A  separate  program  is  provided  to  run  one  turbule  by  itself  * 

*  or  one  tube  of  turbules,  all  the  same  size.  It  is  called  ' slowone . f or ' .  * 

*  Subroutines:  choise,  direct,  dot,  multi,  onetrb,  output,  rotate,  tau.  * 

******************************************************************************* 
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PROGRAM  MAIN 


SUBROUTINE  CHOICE (CONFIG) 

c  This  subroutine  sets  parameters  for  predefined  configurations 
c  chosen  by  the  user  in  the  'input.txt'  file.  It  is  assumed  that 
c  configuration  a  has  been  run  before  configuration  b,  etc. 

SUBROUTINE  DIRECT (AMPTT , AM PVT ) 

c  This  subroutine  calculates  'AMPTT'  and  ' AMPVT '  one  at  a  time  and 
c  accumulates  the  results  directly  (i.e.,  no  time-shifting), 
c  The  'MSIMSL'  below  stands  for  'Microsoft  International  Mathematical 
c  and  Statistical  Library' .  MSIMSL  contains  a  random  number  generator 
c  that  the  'DIRECT'  subroutine  uses. 

SUBROUTINE  DOT (VI , V2 , V1V2 ) 

c  This  subroutine  calculates  the  dot  product  of  two,  real,  three 
c  component  vectors  'VI'  and  'V2',  and  names  this  dot  product  'V1V2' . 

SUBROUTINE  MULTI (AMPT , AM PV) 

c  The  'MULTI'  subroutine  is  a  key  subroutine  in  this  program.  It  is  named 
c  'MULTI'  because  it  is  what  drives  the  calculations  for  multiple  tubes 
c  with  multiple  turbules. 

SUBROUTINE  ONETRB (BO , AT , BV, K) 

c  The  'ONETRB'  subroutine  calculates  'AT',  BV' ,  and  'K' 
c  for  one  turbule,  as  the  turbule  travels  the  total  length  of 
c  the  tube.  The  'ONETRB'  subroutine  is  used  by  both  the  'DIRECT' 
c  subroutine,  and  the  time -shifting  'TAU'  subroutine. 

c  The  'MSIMSL'  subroutine  library  is  need  for  the  random  number  generator. 
SUBROUTINE  OUTPUT (AMPT , AMPV) 

c  This  subroutine  takes  the  complex  amplitudes  'AMPT'  and  'AMPV' 
c  as  functions  of  time,  and  calculates  the  magnitude,  and  the  real 
c  and  imaginary  parts.  It  also  performs  the  FFT  on  them.  Then  the 
c  subroutine  writes  the  data  to  files  called  'time.txt'  and  'fft.txt' 

SUBROUTINE  ROTATE (ALPHA, MW) 

c  The  'ROTATE'  subroutine  rotates  the  coordinate  system  so  that  the  wind 
c  direction  is  along  the  x-axis.  The  origin  is  also  moved  to  the  point 
c  located  half  way  between  the  source  and  detector. 

SUBROUTINE  TAU (AMPTT , AMPVT) 

c  The  'TAU'  subroutine  employs  the  time  -  shifting  procedure. 


* 


* 


*  -  List  of  Symbols 

*  _ 


*  A 

*  AAT(NMAX) 

* 

*  AAV  (NMAX) 

* 

*  AF 

*  ALPHA 

*  AMPT(NMAX) 

*  AMPTT (NMAX) 

* 

*  AMPV (NMAX) 

*  AMPVT (NMAX) 


real;  characteristic  size  of  turbule. 

complex;  temperature  amplitude  accumulation  variable  used 
in  the  'CHOICE'  subroutine. 

complex;  velocity  amplitude  accumulation  variable  used 
in  the  'CHOICE'  subroutine. 

real;  'angular  frequency'  equal  to  2*Pi*nu_s. 
real;  angle  of  wind  relative  to  non-rotated  x-axis. 

complex;  amplitude  due  to  temperature  fluctuations, 
complex;  temperature  amplitude  accumulation  variable  for  a 
tube  used  in  the  'DIRECT'  and  'TAU'  subroutines, 
complex;  amplitude  due  to  velocity  fluctuations, 
complex;  velocity  amplitude  accumulation  variable  for  a 


★ 

★ 
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*  AT(NMAX)  : 

★ 

*  ATO  (NMAX)  : 

* 

*  AV  (NMAX)  : 

* 

*  AVO  : 

* 

*  AP  : 

*  A1  : 

*  BCY  : 

*  BCZ  : 

*  BRK  : 

*  BT  : 

*  BTA ( 3 )  : 

* 

* 

*  BTARS  : 

*  BTARO  : 

*  BV { NMAX , 3 )  : 

*  BO (3 ) 

*  C  : 

*  CALPHA  : 

*  COEF  : 

*  CONFIG  : 

* 

*  CPHI  : 

*  CROSS (3)  : 

*  CTEMP  : 

*  C THETA  : 

*  D  (3  )  : 

*  DA  : 

*  DCM (2 )  : 

* 

*  DI 

*  DNS  : 

*  DS  : 

*  DT  : 

*  DTMAX  : 

*  DW ( 2 ) 

* 

*  Dl  : 

*  ETA  (3)  : 

*  FAMPT (NMAX) : 

*  FAMPV (NMAX) : 

*  FRQ(NMAX)  : 

* 

*  GAMMA2  : 

* 

*  GAMT  : 

*  GAMV  : 

*  GAUSS  : 

* 

*  I  : 

*  IAMPT  : 

*  IAMPV  : 

*  II  : 

*  IM  : 


tube  used  in  the  'DIRECT'  and  'TAU'  subroutines, 
complex;  temperature  amplitude  for  one  turbule  used 
in  the  ' ONETRB '  subroutine. 

complex;  temperature  amplitude  for  one  turbule  used 
in  the  'TAU'  subroutine  in  the  far  field, 
complex;  velocity  amplitude  for  one  turbule  used 
in  the  'ONETRB'  subroutine. 

complex;  velocity  amplitude  for  one  turbule  used 
in  the  'TAU'  subroutine  in  the  far  field, 
real;  'a  prime' ,  reference  distance  from  the  source, 
real;  size  of  largest  turbule. 

real;  y  coordinate  of  centerline  of  largest  tube, 
real;  z  coordinate  of  centerline  of  largest  tube, 
complex;  'bracket',  a  defined  variable. 

real;  the  location  of  the  turbule  as  a  function  of  time. 

real;  the  Greek  letter  'beta'.  One  of  the 

IMSL  library  subroutines  used  'beta',  so  'BTA'  was  used 

as  an  alternate  abbreviation. 

real;  the  dot  product  of  'BTA'  and  'RS'. 

real;  the  dot  product  of  'BTA'  and  'RO' . 

complex;  'B  sub  v' ,  a  defined  variable. 

real;  initial  location  of  turbule. 

real;  speed  of  sound. 

real;  'cosine  of  alpha'. 

real;  'coefficient',  a  defined  variable. 

character;  abbreviation  for  'configuration'  used  to  choose 
predefined  tube  configurations, 
real;  'cosine  of  phi'. 

real;  the  cross  product  of  ' RDTH '  and  'K'. 

complex;  a  temporary  complex  variable. 

real;  'cosine  of  theta'. 

real;  position  vector  to  the  detector. 

real;  distance  between  center  of  cells  for  size  A. 

real;  position  vector  to  detector  relative  to  'center  of  mass' 

position  of  source  and  detector. 

integer;  time -shifting  increment  index. 

real;  distance  between  center  of  cells  for  smallest  size. 

real;  'D  sub  s',  a  defined  variable. 

real;  time  step  increment. 

real;  maximum  time  step  increment. 

real;  x  and  y  components  of  position  vector  to  detector 
in  the  rotated  'wind'  frame. 

real;  distance  between  center  of  cells  for  largest  size  A1 . 
real;  Greek  letter  'eta',  random  displacement  vector, 
complex;  Fourier  transform  of  'AMPT' . 
complex;  Fourier  transform  of  'AMPV' . 

real;  abbreviation  for  'frequency',  the  independent  variable 
of  the  Fourier  transform. 

real;  abbreviation  for  'lower  case  gamma,  squared', 

a  defined  variable. 

real;  'capital  gamma  sub  T' . 

real;  'capital  gamma  sub  v' . 

real;  'gaussian',  a  variable  that  has  a  form  that 
resembles  a  gaussian  function, 
integer;  time  index,  with  range  from  1  to  N. 
real;  'imaginary  part  of  AMPT' . 
real;  'imaginary  part  of  AMPV' . 

complex;  the  imaginary  number  'i',  the  square  root  of  -1. 
integer;  index  for  turbule  in  a  tube,  with  range  from  1  to  M. 
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*  INS 

*  IS 

*  IT 

*  II 

*  J 

*  J2 

*  J3 

*  K ( NMAX , 3 ) 

*  KETA 

*  KS 

*  L 

*  LMIN 

*  M 

*  MAMPT 

*  MAMPV 

*  MBTA 

*  MFAMPT 

*  MFAMPV 

*  MK 

*  MRDT 

*  MRS 

*  MRTS 

*  MRO 

*  MW 

*  Ml 

*  N 

*  NMAX 

* 

*  NMIN 

★ 

*  NS 

*  NST 

*  NT 

*  NUS 

* 

*  OMEGAH ( 3 ) 

* 

*  OMEGA1 

*  PF 

*  PHASE 

*  pj 

*  RAMPT 

*  RAMPV 

*  RDT { 3 ) 

*  RDTH  ( 3 ) 

*  RN  (4  ) 

*  RNA  {  3  ) 

*  RNT  ( 2  ) 

*  RPHASE 

*  RS  ( 3  ) 

*  RTS (3) 

*  RTSH ( 3 ) 

* 

*  RO (3) 

*  ROH  {  3  ) 

* 

*  RORDTH 

*  S  (3) 

*  SALPHA 


:  integer;  number  of  time  points  in  smallest  cell. 

:  integer;  index  for  sizes  with  range  from  1  to  NS . 

:  integer;  index  for  largest  tube,  with  range  from  1  to  NT. 

:  integer;  number  of  time  points  in  largest  cell. 

:  integer;  index  for  vector  components,  with  range  1  to  3 . 

:  integer;  running  index  for  y  component  of  a  vector. 

:  integer;  running  index  for  z  component  of  a  vector. 

:  real;  a  defined  variable. 

:  real;  the  dot  product  of  '  K'  and  'ETA' . 

:  real;  wave  number  for  the  source  'k  sub  s' . 

:  real;  length  of  tubes. 

:  real;  minimum  length  that  tubes  can  have. 

:  integer;  number  of  turbules  in  a  tube. 

:  real;  'magnitude  of  AMPT' . 

:  real;  'magnitude  of  AMPV' . 

:  real;  'magnitude  of  BTA' . 

:  real;  'magnitude  of  Fourier  transform  of  AMPT' . 

:  real;  'magnitude  of  Fourier  transform  of  AMPV'. 

:  real;  'magnitude  of  K' . 

:  real;  'magnitude  of  RDT' . 

:  real;  'magnitude  of  RS' . 

:  real;  'magnitude  of  RTS'  . 

:  real;  'magnitude  of  RO'. 

:  real;  'magnitude  of  W' . 

:  integer;  number  of  turbules  in  largest  tube. 

:  integer;  number  of  time  (and  frequency)  data  points  used. 

:  integer;  maximum  number  of  time  (and  frequency)  data  points 
that  can  be  used. 

:  integer;  minimum  number  of  time  (and  frequency)  data  points 
that  can  be  used. 

:  integer;  number  of  sizes  of  turbules. 

:  integer;  number  of  tubes  in  largest  cell  after  each  halving. 

:  integer;  number  of  largest  tubes. 

:  real;  'Greek  letter  nu  sub  s',  the  frequency 
of  the  source . 

:  real;  the  unit  vector  'omega  hat'.  Omega 
is  rotation  velocity  vector  of  a  turbule. 

:  real;  scaling  parameter  for  the  velocity  turbules. 

:  real;  'packing  fraction'. 

:  real;  a  defined  variable  in  a  complex  exponent. 

:  real;  variable  for  3.14159... 

:  real;  'real  part  of  AMPT' . 

:  real;  'real  part  of  AMPV' . 

:  real;  vector  from  turbule  to  detector. 

:  real;  'RDT  hat',  unit  vector  in  direction  of  'RDT'. 

:  real;  a  random  number  array. 

:  real;  a  random  number  array  for  angles. 

:  real;  a  random  number  array  for  temperature  sign  of  +  or  - . 

:  complex;  abbreviation  for  the  random  phase,  a  defined  variable. 

:  real;  a  defined  variable. 

:  real;  vector  from  source  to  turbule. 

:  real;  abbreviation  for  'RTS  hat',  a  unit  vector  in  the 
direction  of  'RTS' . 

:  real;  a  defined  variable. 

:  real;  abbreviation  for  ' RO  hat',  a  unit  vector  in  the 
direction  of  ' RO  '  . 

:  real;  the  dot  product  of  'ROH'  and  'RDTH' . 

:  real;  position  vector  to  the  source. 

:  real;  'sine  of  alpha'. 
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*  SCM  (2 ) 

* 

*  SIGMA2 

*  SPHI 

*  S THETA 

*  SW(2) 

* 

*  T 

*  TEMP 

*  TIMER 

*  TMAX 

*  TS 

*  TSF 

*  TWOPI 

*  VSF 

*  VI  ( 3  ) 

*  V2 ( 3 } 

*  V1V2 

* 

*  W  {  3  ) 

*  XDS 

* 

*  YC 

*  YC1 (20) 

*  ZB 

*  zc 

*  ZC1 (20) 

*  ZI(3) 


:  real;  position  vector  to  source  relative  to  'center  of  mass' 
position  of  source  and  detector. 

:  real;  'sigma  squared',  a  defined  variable. 

:  real;  'sine  of  phi'. 

:  real;  'sine  of  theta' . 

:  real;  x  and  y  components  of  position  vector  to  source 
in  the  rotated  'wind'  frame. 

:  real;  time  variable. 

:  real;  a  temporary  real  variable. 

:  real;  gives  value  of  total  time  for  the  program  to  run. 

:  real;  largest  value  of  time  variable. 

:  real;  time  for  calculations  for  one  size  to  be  completed. 

:  real;  'temperature  scale  factor',  a  defined  variable. 

;  real ;  ' 2*pi ' . 

:  real;  'velocity  scale  factor',  a  defined  variable. 

:  real;  arbitrary  vector  used  in  the  'DOT'  subroutine. 

:  real;  arbitrary  vector  used  in  the  'DOT'  subroutine. 

:  real;  result  of  the  dot  product  of  'VI'  and  'V2'  in  the 

'DOT'  subroutine 
:  real;  wind  velocity  vector. 

:  real;  x  component  of  the  vector  from  the  source 
to  the  detector. 

:  real;  y  coordinate  to  center  line  of  tube  for  size  A. 

:  real;  y  coordinate  to  center  line  of  tube  for  size  A1 . 

:  real;  height  above  ground  where  scattering  volume  begins. 

:  real;  z  coordinate  to  center  line  of  tube  for  size  A. 

:  real;  z  coordinate  to  center  line  of  tube  for  size  A1 . 

:  real;  Greek  letter  'zi',  random  variable 
between  +1/4  and  -1/4. 


Descriptions  of  Program  Inputl  and  Program  Slow  One  follow: 

******************************************************************************* 

*  Program  Inputl  * 

*  Version  1,0:  11-21-98  * 

**********  *****************************************************  **************** 

*  This  program  enters  the  environment  parameters,  and  the  turbule  * 

*  parameters  needed  for  the  slowone . f or  program.  These  parameters  * 

*  are  written  to  the  file  " inputl . txt . "  * 
******************************************************************************* 


******************  **************  ********  ********************************  ******* 

*  Program  Slow  One  * 

*  Version  1.0:  11-21-98  * 

******************************************************************************* 

*  This  program  calculates  the  pressure  amplitudes  as  a  function  of  time,  * 

*  and  their  Fourier  transforms  as  a  function  of  frequency.  This  program  * 

*  uses  input  from  the  file  'inputl.txt'  that  is  generated  by  running  the  * 

*  'inputl. for'  Fortran  program.  The  Slow  One  program  is  designed  to  run  * 

*  one  tube  of  turbules  of  one  size,  or  just  one  turbule  by  itself.  * 

*  Subroutines:  choice,  direct,  dot,  multi,  onetrb,  output,  rotate,  tau.  * 

********************************************************************************** 
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Subroutines  and  symbols  for  Program  Inputl /Program  Slow  One  are  al¬ 
most  identical  with  those  of  Program  Input/Program  Slow  Eddy. 

Flow  charts  for  Program  Slow  Eddy  and  its  eight  subroutines  are  given  in 
figures  D-l,  D-2,  D-3,  D-4,  and  D-5. 


Figure  D-l.  Flow  chart 
for  main  program. 


/  Read 
/  Input 
/  File 
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Figure  D-2.  Flow  chart 
for  subroutine  CHOICE. 
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Figure  D-5.  Flow  chart 
for  subroutines 
ONETRB,  OUTPUT, 
and  DOT. 


Not  shown  in  figure  D-2  for  subroutine  CHOICE  is  the  intermediate  save 
operations  in  branches  for  which  the  variable  CONFIG  is  b  through  e.  These 
are  different  in  each  branch,  with  typical  file  names  given  in  table  D-l. 


Table  D-l.  Intermediate 
save  operations. 


CONFIG 

File  name(s) 

b 

cfl2.txt 

c 

cf34.txt  and  cfl234.txt 

d 

cf56.txt 

e 

cfl23456.txt* 

*This  branch  does  not  call  subroutine 
MULTI. 
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Distribution 
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The  Pentagon 

Washington  DC  20301-7100 

Ofc  of  the  Secy  of  Defns 
Attn  OUSD(A&T)/ODDR&E(R)  RJ  Trew 
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